From 81dbff308753b8e0ee543959ab77e364b686bd59 Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Wed, 15 Jan 2025 16:15:50 +0100
Subject: [PATCH] Move everything into the Dune::GFE namespace

Previously, some code was in Dune::GFE, some was in Dune::,
and a lot of it was in the global namespace.
---
 CHANGELOG.md                                  |  2 +
 dune/gfe/assemblers/cosseratrodenergy.hh      |  3 +-
 dune/gfe/assemblers/geodesicfeassembler.hh    |  6 +-
 .../assemblers/geodesicfeassemblerwrapper.hh  |  6 +-
 .../gfe/assemblers/l2distancesquaredenergy.hh |  6 ++
 dune/gfe/assemblers/localenergy.hh            | 13 ++--
 dune/gfe/assemblers/localfirstordermodel.hh   |  2 +-
 .../localgeodesicfeadolcstiffness.hh          |  6 ++
 .../assemblers/localgeodesicfefdstiffness.hh  |  6 ++
 .../assemblers/localgeodesicfestiffness.hh    | 13 ++--
 dune/gfe/assemblers/localintegralenergy.hh    |  3 +-
 dune/gfe/assemblers/mixedgfeassembler.hh      |  5 ++
 .../nonplanarcosseratshellenergy.hh           |  5 +-
 dune/gfe/assemblers/simofoxenergy.hh          |  5 +-
 dune/gfe/assemblers/sumenergy.hh              |  4 +-
 .../surfacecosseratstressassembler.hh         |  8 ++-
 dune/gfe/assemblers/weightedsumenergy.hh      |  6 ++
 dune/gfe/averagedistanceassembler.hh          | 10 ++-
 dune/gfe/cosseratstrain.hh                    |  9 +--
 dune/gfe/cosseratvtkreader.hh                 |  8 +--
 dune/gfe/cosseratvtkwriter.hh                 |  5 ++
 dune/gfe/densities/bulkcosseratdensity.hh     |  3 +-
 dune/gfe/densities/localdensity.hh            |  3 +-
 .../densities/planarcosseratshelldensity.hh   |  2 +-
 dune/gfe/filereader.hh                        |  7 +-
 .../functions/embeddedglobalgfefunction.hh    |  4 +-
 dune/gfe/functions/globalgfefunction.hh       |  4 +-
 dune/gfe/functions/localgeodesicfefunction.hh | 19 +++--
 .../gfe/functions/localprojectedfefunction.hh |  8 +--
 .../functions/localquickanddirtyfefunction.hh |  8 +--
 dune/gfe/geodesicdifference.hh                |  6 ++
 dune/gfe/geodesicfefunctionadaptor.hh         |  6 ++
 dune/gfe/geometries/moebiusstrip.hh           |  5 +-
 dune/gfe/geometries/twistedstrip.hh           |  5 +-
 dune/gfe/gfedifferencenormsquared.hh          |  6 ++
 dune/gfe/gramschmidtsolver.hh                 | 14 ++--
 dune/gfe/linearalgebra.hh                     |  8 +--
 dune/gfe/maxnormtrustregion.hh                |  6 ++
 dune/gfe/mixedriemanniantrsolver.cc           |  4 +-
 dune/gfe/mixedriemanniantrsolver.hh           |  6 ++
 dune/gfe/neumannenergy.hh                     |  4 +-
 dune/gfe/parallel/globalmapper.hh             |  6 +-
 dune/gfe/parallel/globalp1mapper.hh           |  6 +-
 dune/gfe/parallel/globalp2mapper.hh           |  9 ++-
 dune/gfe/parallel/mapperfactory.hh            |  4 +-
 dune/gfe/parallel/matrixcommunicator.hh       |  5 ++
 dune/gfe/parallel/mpifunctions.hh             |  5 ++
 dune/gfe/parallel/p2mapper.hh                 |  6 ++
 dune/gfe/parallel/vectorcommunicator.hh       |  5 ++
 dune/gfe/polardecomposition.hh                |  6 +-
 dune/gfe/riemannianpnsolver.cc                | 10 +--
 dune/gfe/riemannianpnsolver.hh                |  9 ++-
 dune/gfe/riemanniantrsolver.cc                | 22 +++---
 dune/gfe/riemanniantrsolver.hh                |  9 ++-
 dune/gfe/rodfactory.hh                        |  6 ++
 dune/gfe/skewmatrix.hh                        |  5 ++
 dune/gfe/spaces/hyperbolichalfspacepoint.hh   |  6 ++
 dune/gfe/spaces/orthogonalmatrix.hh           |  6 ++
 dune/gfe/spaces/productmanifold.hh            |  2 +-
 dune/gfe/spaces/quaternion.hh                 |  8 ++-
 dune/gfe/spaces/realtuple.hh                  |  9 ++-
 dune/gfe/spaces/rotation.hh                   | 12 ++--
 dune/gfe/spaces/unitvector.hh                 | 12 +++-
 dune/gfe/svd.hh                               |  5 ++
 dune/gfe/symmetricmatrix.hh                   |  6 +-
 dune/gfe/targetspacertrsolver.cc              |  4 +-
 dune/gfe/targetspacertrsolver.hh              |  8 ++-
 dune/gfe/tensor3.hh                           |  6 +-
 dune/gfe/tensorssd.hh                         |  5 ++
 dune/gfe/trustregionmmgbasesolver.hh          |  4 ++
 dune/gfe/vtkfile.hh                           |  9 +--
 src/compute-disc-error.cc                     | 72 +++++++++----------
 src/cosserat-continuum.cc                     | 58 +++++++--------
 src/cosserat-rod.cc                           | 18 ++---
 src/film-on-substrate-stress-plot.cc          |  8 +--
 src/film-on-substrate.cc                      | 40 +++++------
 src/gradient-flow.cc                          | 27 +++----
 src/harmonicmaps.cc                           | 25 +++----
 src/simofoxshell.cc                           | 18 ++---
 test/averagedistanceassemblertest.cc          | 18 ++---
 test/cosseratcontinuumtest.cc                 | 26 +++----
 test/cosseratrodenergytest.cc                 | 14 ++--
 test/cosseratrodtest.cc                       | 14 ++--
 test/embeddedglobalgfefunctiontest.cc         |  4 +-
 test/filmonsubstratetest.cc                   | 42 +++++------
 test/harmonicmaptest.cc                       |  8 +--
 test/interpolationderivativestest.cc          | 12 ++--
 test/localadolcstiffnesstest.cc               | 12 ++--
 test/localgeodesicfefunctiontest.cc           | 52 +++++++-------
 test/localintegralstiffnesstest.cc            | 64 ++++++++---------
 test/localprojectedfefunctiontest.cc          | 44 ++++++------
 test/mixedriemannianpnsolvertest.cc           | 26 +++----
 test/multiindex.hh                            |  6 ++
 test/nonplanarcosseratenergytest.cc           | 10 +--
 test/orthogonalmatrixtest.cc                  | 12 ++--
 test/planarcosseratshelltest.cc               | 30 ++++----
 test/polardecompositiontest.cc                |  2 +-
 test/rotationtest.cc                          | 56 +++++++--------
 test/surfacecosseratstressassemblertest.cc    | 22 +++---
 test/svdtest.cc                               |  4 +-
 test/targetspacetest.cc                       | 28 ++++----
 test/valuefactory.hh                          |  5 ++
 102 files changed, 731 insertions(+), 520 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 9d932555..378c61e9 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,7 @@
 # Master
 
+- The entire code is now in the `Dune::GFE` namespace.
+
 - Remove the files `localgfetestfunctionbasis.hh` and `localtangentfefunction.hh`.
   They were never used for anything at all, I and currently wouldn't know
   what to use them for, either.
diff --git a/dune/gfe/assemblers/cosseratrodenergy.hh b/dune/gfe/assemblers/cosseratrodenergy.hh
index 2cb9fb43..7d6c6e67 100644
--- a/dune/gfe/assemblers/cosseratrodenergy.hh
+++ b/dune/gfe/assemblers/cosseratrodenergy.hh
@@ -18,7 +18,8 @@
 #include <dune/gfe/spaces/realtuple.hh>
 #include <dune/gfe/spaces/rotation.hh>
 
-namespace Dune::GFE {
+namespace Dune::GFE
+{
 
   template<class Basis, class LocalInterpolationRule, class RT>
   class CosseratRodEnergy
diff --git a/dune/gfe/assemblers/geodesicfeassembler.hh b/dune/gfe/assemblers/geodesicfeassembler.hh
index b24be867..0d976dc1 100644
--- a/dune/gfe/assemblers/geodesicfeassembler.hh
+++ b/dune/gfe/assemblers/geodesicfeassembler.hh
@@ -10,6 +10,10 @@
 
 #include <dune/solvers/common/wrapownshare.hh>
 
+
+namespace Dune::GFE
+{
+
 /** \brief A global FE assembler for problems involving functions that map into non-Euclidean spaces
  */
 template <class Basis, class TargetSpace>
@@ -273,6 +277,6 @@ computeEnergy(const std::vector<TargetSpace>& sol) const
 
 }
 
-
+}  // namespace Dune::GFE
 
 #endif
diff --git a/dune/gfe/assemblers/geodesicfeassemblerwrapper.hh b/dune/gfe/assemblers/geodesicfeassemblerwrapper.hh
index c751263c..692e963c 100644
--- a/dune/gfe/assemblers/geodesicfeassemblerwrapper.hh
+++ b/dune/gfe/assemblers/geodesicfeassemblerwrapper.hh
@@ -5,7 +5,8 @@
 #include <dune/common/tuplevector.hh>
 #include <dune/gfe/spaces/productmanifold.hh>
 
-namespace Dune::GFE {
+namespace Dune::GFE
+{
 
   /** \brief A wrapper that wraps a MixedGFEAssembler into an assembler that does not distinguish between the two finite element spaces
 
@@ -72,7 +73,8 @@ namespace Dune::GFE {
     auto splitVector(const std::vector<TargetSpace>& sol) const;
     std::unique_ptr<MatrixType> hessianMixed_;
   }; // end class
-} //end namespace
+
+}  // namespace Dune::GFE
 
 
 template <class Basis, class ScalarBasis, class TargetSpace>
diff --git a/dune/gfe/assemblers/l2distancesquaredenergy.hh b/dune/gfe/assemblers/l2distancesquaredenergy.hh
index d2b55d4c..d5b5a5c1 100644
--- a/dune/gfe/assemblers/l2distancesquaredenergy.hh
+++ b/dune/gfe/assemblers/l2distancesquaredenergy.hh
@@ -9,6 +9,10 @@
 #include <dune/gfe/functions/globalgfefunction.hh>
 #include <dune/gfe/functions/localgeodesicfefunction.hh>
 
+
+namespace Dune::GFE
+{
+
 template<class Basis, class TargetSpace>
 class L2DistanceSquaredEnergy
   : public Dune::GFE::LocalEnergy<Basis,TargetSpace>
@@ -87,4 +91,6 @@ public:
   }
 };
 
+}  // namespace Dune::GFE
+
 #endif
diff --git a/dune/gfe/assemblers/localenergy.hh b/dune/gfe/assemblers/localenergy.hh
index 9056e5b7..57cbee35 100644
--- a/dune/gfe/assemblers/localenergy.hh
+++ b/dune/gfe/assemblers/localenergy.hh
@@ -7,9 +7,9 @@
 
 #include <dune/gfe/spaces/productmanifold.hh>
 
-namespace Dune {
-
-  namespace GFE:: Impl
+namespace Dune::GFE
+{
+  namespace Impl
   {
     /** \brief A class exporting container types for coefficient sets
      *
@@ -34,9 +34,8 @@ namespace Dune {
       using Coefficients = std::vector<ProductManifold<Factors...> >;
       using CompositeCoefficients = TupleVector<std::vector<Factors>... >;
     };
-  }
+  }  // namespace Impl
 
-  namespace GFE {
 
     /** \brief Base class for energies defined by integrating over one grid element */
     template<class Basis, class TargetSpace>
@@ -69,8 +68,6 @@ namespace Dune {
 
     };
 
-  } // namespace GFE
-
-}  // namespace Dune
+}  // namespace Dune::GFE
 
 #endif  // DUNE_GFE_LOCALENERGY_HH
diff --git a/dune/gfe/assemblers/localfirstordermodel.hh b/dune/gfe/assemblers/localfirstordermodel.hh
index 7ab06435..30fc7ca7 100644
--- a/dune/gfe/assemblers/localfirstordermodel.hh
+++ b/dune/gfe/assemblers/localfirstordermodel.hh
@@ -20,6 +20,6 @@ namespace Dune::GFE
 
   };
 
-}  // namespace Dune
+}  // namespace Dune::GFE
 
 #endif   // DUNE_GFE_LOCALFIRSTORDERMODEL_HH
diff --git a/dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh b/dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh
index 7b16b831..cda65b84 100644
--- a/dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh
+++ b/dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh
@@ -14,6 +14,10 @@
 #include <dune/gfe/assemblers/localgeodesicfestiffness.hh>
 #include <dune/gfe/spaces/productmanifold.hh>
 
+
+namespace Dune::GFE
+{
+
 /** \brief Assembles energy gradient and Hessian with ADOL-C (automatic differentiation)
  */
 template<class Basis, class TargetSpace>
@@ -825,4 +829,6 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
   }
 }
 
+}  // namespace Dune::GFE
+
 #endif
diff --git a/dune/gfe/assemblers/localgeodesicfefdstiffness.hh b/dune/gfe/assemblers/localgeodesicfefdstiffness.hh
index 630e9a1e..e08c76af 100644
--- a/dune/gfe/assemblers/localgeodesicfefdstiffness.hh
+++ b/dune/gfe/assemblers/localgeodesicfefdstiffness.hh
@@ -6,6 +6,10 @@
 
 #include <dune/gfe/assemblers/localgeodesicfestiffness.hh>
 
+
+namespace Dune::GFE
+{
+
 /** \brief Assembles energy gradient and Hessian with finite difference approximations
  */
 template<class Basis, class TargetSpace, class field_type=double>
@@ -274,4 +278,6 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
   }
 }
 
+}  // namespace Dune::GFE
+
 #endif
diff --git a/dune/gfe/assemblers/localgeodesicfestiffness.hh b/dune/gfe/assemblers/localgeodesicfestiffness.hh
index e1f936fa..b74769b6 100644
--- a/dune/gfe/assemblers/localgeodesicfestiffness.hh
+++ b/dune/gfe/assemblers/localgeodesicfestiffness.hh
@@ -64,8 +64,7 @@ namespace Dune::GFE
       // Type of the local Hessian
       using CompositeHessian = FieldMatrix<Matrix<double>,2,2>;
     };
-  }
-}
+  }  // namespace Impl
 
 template<class Basis, class TargetSpace>
 class LocalGeodesicFEStiffness
@@ -77,16 +76,18 @@ public:
 
    */
   virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
-                                          const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Coefficients& coefficients,
+                                          const typename GFE::Impl::LocalStiffnessTypes<TargetSpace>::Coefficients& coefficients,
                                           std::vector<double>& localGradient,
-                                          typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Hessian& localHessian) const = 0;
+                                          typename GFE::Impl::LocalStiffnessTypes<TargetSpace>::Hessian& localHessian) const = 0;
 
   /** \brief Assemble the local gradient and stiffness matrix at the current position -- Composite version
    */
   virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
-                                          const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeCoefficients& coefficients,
+                                          const typename GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeCoefficients& coefficients,
                                           std::vector<double>& localGradient,
-                                          typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeHessian& localHessian) const = 0;
+                                          typename GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeHessian& localHessian) const = 0;
 };
 
+}  // namespace Dune::GFE
+
 #endif
diff --git a/dune/gfe/assemblers/localintegralenergy.hh b/dune/gfe/assemblers/localintegralenergy.hh
index 4424a731..8458ad0e 100644
--- a/dune/gfe/assemblers/localintegralenergy.hh
+++ b/dune/gfe/assemblers/localintegralenergy.hh
@@ -10,7 +10,8 @@
 #include <dune/gfe/densities/localdensity.hh>
 
 
-namespace Dune::GFE {
+namespace Dune::GFE
+{
 
 #if ! DUNE_VERSION_GTE(DUNE_LOCALFUNCTIONS, 2, 10)
   namespace Impl
diff --git a/dune/gfe/assemblers/mixedgfeassembler.hh b/dune/gfe/assemblers/mixedgfeassembler.hh
index 67c484ab..e68d5a08 100644
--- a/dune/gfe/assemblers/mixedgfeassembler.hh
+++ b/dune/gfe/assemblers/mixedgfeassembler.hh
@@ -12,6 +12,9 @@
 #include <dune/gfe/assemblers/localgeodesicfestiffness.hh>
 
 
+namespace Dune::GFE
+{
+
 /** \brief A global FE assembler for problems involving functions that map into non-Euclidean spaces
  */
 template <class Basis, class TargetSpace>
@@ -341,4 +344,6 @@ computeEnergy(const std::vector<TargetSpace0>& configuration0,
 
 }
 
+}  // namespace Dune::GFE
+
 #endif
diff --git a/dune/gfe/assemblers/nonplanarcosseratshellenergy.hh b/dune/gfe/assemblers/nonplanarcosseratshellenergy.hh
index f2110a6f..15e55ceb 100644
--- a/dune/gfe/assemblers/nonplanarcosseratshellenergy.hh
+++ b/dune/gfe/assemblers/nonplanarcosseratshellenergy.hh
@@ -60,8 +60,7 @@ namespace Dune::GFE
         return localView.tree().finiteElement();
       }
     };
-  }
-}
+  }  // namespace Impl
 
 /** \brief Assembles the cosserat energy for a single element.
  *
@@ -348,4 +347,6 @@ energy(const typename Basis::LocalView& localView,
   return energy;
 }
 
+}  // namespace Dune::GFE
+
 #endif   //#ifndef DUNE_GFE_NONPLANARCOSSERATSHELLENERGY_HH
diff --git a/dune/gfe/assemblers/simofoxenergy.hh b/dune/gfe/assemblers/simofoxenergy.hh
index 3dd02a96..a59fc9e4 100644
--- a/dune/gfe/assemblers/simofoxenergy.hh
+++ b/dune/gfe/assemblers/simofoxenergy.hh
@@ -12,7 +12,8 @@
 #include <dune/gfe/spaces/realtuple.hh>
 #include <dune/gfe/spaces/unitvector.hh>
 
-namespace Dune::GFE {
+namespace Dune::GFE
+{
 
   /** \brief Get LocalFiniteElements from a localView, for different tree depths of the local view
    *
@@ -360,5 +361,7 @@ namespace Dune::GFE {
     }
     return energy;
   }
+
 }  // namespace Dune::GFE
+
 #endif  // DUNE_GFE_SIMOFOX_ENERGY_HH
diff --git a/dune/gfe/assemblers/sumenergy.hh b/dune/gfe/assemblers/sumenergy.hh
index 2488d273..3f8a43ac 100644
--- a/dune/gfe/assemblers/sumenergy.hh
+++ b/dune/gfe/assemblers/sumenergy.hh
@@ -5,8 +5,8 @@
 
 #include <dune/gfe/assemblers/localenergy.hh>
 
-namespace Dune::GFE {
-
+namespace Dune::GFE
+{
   /**
      \brief Assembles the a sum of energies for a single element by summing up the energies of each GFE::LocalEnergy.
 
diff --git a/dune/gfe/assemblers/surfacecosseratstressassembler.hh b/dune/gfe/assemblers/surfacecosseratstressassembler.hh
index bc451e24..557488d7 100644
--- a/dune/gfe/assemblers/surfacecosseratstressassembler.hh
+++ b/dune/gfe/assemblers/surfacecosseratstressassembler.hh
@@ -9,7 +9,9 @@
 
 #include <dune/matrix-vector/transpose.hh>
 
-namespace Dune::GFE {
+namespace Dune::GFE
+{
+
   /** \brief An assembler that can calculate the norms of specific stress tensors for each element for an output by film-on-substrate
 
      \tparam BasisOrderD Basis used for the displacement
@@ -301,5 +303,7 @@ namespace Dune::GFE {
       }
     }
   };
-}
+
+}  // namespace Dune::GFE
+
 #endif
diff --git a/dune/gfe/assemblers/weightedsumenergy.hh b/dune/gfe/assemblers/weightedsumenergy.hh
index bcc189dc..2535be81 100644
--- a/dune/gfe/assemblers/weightedsumenergy.hh
+++ b/dune/gfe/assemblers/weightedsumenergy.hh
@@ -8,6 +8,10 @@
 
 #include <dune/gfe/assemblers/localenergy.hh>
 
+
+namespace Dune::GFE
+{
+
 template<class Basis, class TargetSpace>
 class WeightedSumEnergy
   : public Dune::GFE::LocalEnergy<Basis,TargetSpace>
@@ -54,4 +58,6 @@ public:
   }
 };
 
+}  // namespace Dune::GFE
+
 #endif
diff --git a/dune/gfe/averagedistanceassembler.hh b/dune/gfe/averagedistanceassembler.hh
index 55fd1a33..077483b8 100644
--- a/dune/gfe/averagedistanceassembler.hh
+++ b/dune/gfe/averagedistanceassembler.hh
@@ -5,6 +5,10 @@
 
 #include <dune/gfe/symmetricmatrix.hh>
 
+
+namespace Dune::GFE
+{
+
 /** \tparam TargetSpace The manifold that we are mapping to */
 template <class TargetSpace, class WeightType=double>
 class AverageDistanceAssembler
@@ -70,7 +74,7 @@ public:
   }
 
   void assembleEmbeddedHessian(const TargetSpace& x,
-                               Dune::SymmetricMatrix<ctype,embeddedSize>& matrix) const
+                               SymmetricMatrix<ctype,embeddedSize>& matrix) const
   {
     matrix = 0;
     for (size_t i=0; i<coefficients_.size(); i++)
@@ -81,7 +85,7 @@ public:
   void assembleHessian(const TargetSpace& x,
                        Dune::FieldMatrix<ctype,size,size>& matrix) const
   {
-    Dune::SymmetricMatrix<ctype,embeddedSize> embeddedHessian;
+    SymmetricMatrix<ctype,embeddedSize> embeddedHessian;
     assembleEmbeddedHessian(x,embeddedHessian);
 
     Dune::FieldMatrix<ctype,size,embeddedSize> frame = x.orthonormalFrame();
@@ -98,4 +102,6 @@ public:
 
 };
 
+}  // namespace Dune::GFE
+
 #endif
diff --git a/dune/gfe/cosseratstrain.hh b/dune/gfe/cosseratstrain.hh
index c2405b76..00e96526 100644
--- a/dune/gfe/cosseratstrain.hh
+++ b/dune/gfe/cosseratstrain.hh
@@ -3,9 +3,8 @@
 
 #include <dune/common/fmatrix.hh>
 
-namespace Dune {
-
-  namespace GFE {
+namespace Dune::GFE
+{
 
     /** \brief Strain tensor of a Cosserat material
      *
@@ -107,8 +106,6 @@ namespace Dune {
       FieldMatrix<T,dimworld,dimworld> data_;
     };
 
-  }   // namespace GFE
-
-}  // namespace Dune
+}  // namespace Dune::GFE
 
 #endif   // DUNE_GFE_COSSERATSTRAIN_HH
diff --git a/dune/gfe/cosseratvtkreader.hh b/dune/gfe/cosseratvtkreader.hh
index 659634ff..cda884ee 100644
--- a/dune/gfe/cosseratvtkreader.hh
+++ b/dune/gfe/cosseratvtkreader.hh
@@ -5,10 +5,8 @@
 #include <dune/gfe/spaces/realtuple.hh>
 #include <dune/gfe/spaces/rotation.hh>
 
-namespace Dune
+namespace Dune::GFE
 {
-  namespace GFE
-  {
 
     /** \brief Read configurations of Cosserat models from VTK files into memory */
     class CosseratVTKReader
@@ -39,8 +37,6 @@ namespace Dune
 
     };
 
-  }
-
-}
+}  // namespace Dune::GFE
 
 #endif
diff --git a/dune/gfe/cosseratvtkwriter.hh b/dune/gfe/cosseratvtkwriter.hh
index 69eb2dfb..aaa9f372 100644
--- a/dune/gfe/cosseratvtkwriter.hh
+++ b/dune/gfe/cosseratvtkwriter.hh
@@ -22,6 +22,9 @@
 #include <dune/gfe/spaces/rotation.hh>
 
 
+namespace Dune::GFE
+{
+
 /** \brief Write the configuration of a Cosserat material in VTK format */
 template <class GridView>
 class CosseratVTKWriter
@@ -448,4 +451,6 @@ public:
 
 };
 
+}  // namespace Dune::GFE
+
 #endif
diff --git a/dune/gfe/densities/bulkcosseratdensity.hh b/dune/gfe/densities/bulkcosseratdensity.hh
index 3131cba8..d2a6d975 100644
--- a/dune/gfe/densities/bulkcosseratdensity.hh
+++ b/dune/gfe/densities/bulkcosseratdensity.hh
@@ -10,7 +10,8 @@
 #include <dune/gfe/spaces/realtuple.hh>
 #include <dune/gfe/spaces/rotation.hh>
 
-namespace Dune::GFE {
+namespace Dune::GFE
+{
 
   template<class Element, class field_type>
   class BulkCosseratDensity final
diff --git a/dune/gfe/densities/localdensity.hh b/dune/gfe/densities/localdensity.hh
index f0fd4071..eeb4a5e7 100644
--- a/dune/gfe/densities/localdensity.hh
+++ b/dune/gfe/densities/localdensity.hh
@@ -10,7 +10,8 @@
 
 #include <dune/istl/matrix.hh>
 
-namespace Dune::GFE {
+namespace Dune::GFE
+{
 
   /** \brief A base class for energy densities to be evaluated in an integral energy
    *
diff --git a/dune/gfe/densities/planarcosseratshelldensity.hh b/dune/gfe/densities/planarcosseratshelldensity.hh
index 927448dd..d016626c 100644
--- a/dune/gfe/densities/planarcosseratshelldensity.hh
+++ b/dune/gfe/densities/planarcosseratshelldensity.hh
@@ -296,6 +296,6 @@ namespace Dune::GFE
     double kappa_;
   };
 
-}
+}  // namespace Dune::GFE
 
 #endif   //#ifndef DUNE_GFE_DENSITIES_PLANARCOSSERATSHELLDENSITY_HH
diff --git a/dune/gfe/filereader.hh b/dune/gfe/filereader.hh
index 24efd6ca..88dd87c5 100644
--- a/dune/gfe/filereader.hh
+++ b/dune/gfe/filereader.hh
@@ -8,8 +8,9 @@
 // This is a custom file format parser which reads in a grid deformation file
 // The file must contain lines of the format <grid vertex>:<displacement or rotation>
 // !!! THIS IS A TEMPORARY SOLUTION AND WILL BE REPLACED BY THE Dune::VtkReader in dune-vtk/dune/vtk/vtkreader.hh !!!
-namespace Dune {
-  namespace GFE {
+namespace Dune::GFE
+{
+
     // Convert the pairs {grid vertex, vector of dimension d} in the given file to a map
     template <int d>
     static std::unordered_map<std::string, FieldVector<double,d> > transformFileToMap(std::string pathToFile) {
@@ -38,6 +39,6 @@ namespace Dune {
       }
       return map;
     }
-  }
+
 }
 #endif
diff --git a/dune/gfe/functions/embeddedglobalgfefunction.hh b/dune/gfe/functions/embeddedglobalgfefunction.hh
index 34c821af..531f70b1 100644
--- a/dune/gfe/functions/embeddedglobalgfefunction.hh
+++ b/dune/gfe/functions/embeddedglobalgfefunction.hh
@@ -18,8 +18,8 @@
 
 #include <dune/gfe/functions/globalgfefunction.hh>
 
-namespace Dune::GFE {
-
+namespace Dune::GFE
+{
   template<typename EGGF>
   class EmbeddedGlobalGFEFunctionDerivative;
 
diff --git a/dune/gfe/functions/globalgfefunction.hh b/dune/gfe/functions/globalgfefunction.hh
index cf66c8ed..55357626 100644
--- a/dune/gfe/functions/globalgfefunction.hh
+++ b/dune/gfe/functions/globalgfefunction.hh
@@ -20,8 +20,8 @@
 #include <dune/fufem/functions/virtualgridfunction.hh>
 #endif
 
-namespace Dune::GFE {
-
+namespace Dune::GFE
+{
   namespace Impl {
 
 #if DUNE_VERSION_LTE(DUNE_FUFEM, 2, 9)
diff --git a/dune/gfe/functions/localgeodesicfefunction.hh b/dune/gfe/functions/localgeodesicfefunction.hh
index 6cce0520..0ecd6bb1 100644
--- a/dune/gfe/functions/localgeodesicfefunction.hh
+++ b/dune/gfe/functions/localgeodesicfefunction.hh
@@ -18,6 +18,10 @@
 #include <dune/gfe/tensorssd.hh>
 #include <dune/gfe/linearalgebra.hh>
 
+
+namespace Dune::GFE
+{
+
 // forward declaration
 template <class LocalFiniteElement, class TargetSpace>
 class LocalGfeTestFunctionBasis;
@@ -143,7 +147,7 @@ public:
   }
 private:
 
-  static Dune::SymmetricMatrix<RT,embeddedDim> pseudoInverse(const Dune::SymmetricMatrix<RT,embeddedDim>& dFdq,
+  static SymmetricMatrix<RT,embeddedDim> pseudoInverse(const SymmetricMatrix<RT,embeddedDim>& dFdq,
                                                              const TargetSpace& q)
   {
     const int shortDim = TargetSpace::TangentVector::dimension;
@@ -164,7 +168,7 @@ private:
 
     A.invert();
 
-    Dune::SymmetricMatrix<RT,embeddedDim> result;
+    SymmetricMatrix<RT,embeddedDim> result;
     result = 0.0;
     for (int i=0; i<embeddedDim; i++)
       for (int j=0; j<=i; j++)
@@ -297,7 +301,7 @@ evaluateDerivative(const Dune::FieldVector<ctype, dim>& local, const TargetSpace
 
   AverageDistanceAssembler<TargetSpace> assembler(coefficients_, w);
 
-  Dune::SymmetricMatrix<RT,embeddedDim> dFdq;
+  SymmetricMatrix<RT,embeddedDim> dFdq;
   assembler.assembleEmbeddedHessian(q,dFdq);
 
   // We want to solve
@@ -353,7 +357,7 @@ evaluateDerivativeOfValueWRTCoefficient(const Dune::FieldVector<ctype, dim>& loc
 
   AverageDistanceAssembler<TargetSpace> assembler(coefficients_, w);
 
-  Dune::SymmetricMatrix<RT,embeddedDim> dFdq;
+  SymmetricMatrix<RT,embeddedDim> dFdq;
   assembler.assembleEmbeddedHessian(q,dFdq);
 
   const int shortDim = TargetSpace::TangentVector::dimension;
@@ -461,7 +465,7 @@ evaluateDerivativeOfGradientWRTCoefficient(const Dune::FieldVector<ctype, dim>&
   AverageDistanceAssembler<TargetSpace> assembler(coefficients_, w);
 
   /** \todo Use a symmetric matrix here */
-  Dune::SymmetricMatrix<RT,embeddedDim> dFdq;
+  SymmetricMatrix<RT,embeddedDim> dFdq;
   assembler.assembleEmbeddedHessian(q,dFdq);
 
 
@@ -476,7 +480,7 @@ evaluateDerivativeOfGradientWRTCoefficient(const Dune::FieldVector<ctype, dim>&
   // dFDq is not invertible, if the target space is embedded into a higher-dimensional
   // Euclidean space.  Therefore we use its pseudo inverse.  I don't think that is the
   // best way, though.
-  Dune::SymmetricMatrix<RT,embeddedDim> dFdqPseudoInv = pseudoInverse(dFdq,q);
+  SymmetricMatrix<RT,embeddedDim> dFdqPseudoInv = pseudoInverse(dFdq,q);
 
   //
   Tensor3<RT,embeddedDim,embeddedDim,embeddedDim> dvDqF
@@ -497,7 +501,7 @@ evaluateDerivativeOfGradientWRTCoefficient(const Dune::FieldVector<ctype, dim>&
   TensorSSD<RT, embeddedDim,embeddedDim> dqdwF(coefficients_.size());
 
   for (size_t k=0; k<coefficients_.size(); k++) {
-    Dune::SymmetricMatrix<RT,embeddedDim> hesse = TargetSpace::secondDerivativeOfDistanceSquaredWRTSecondArgument(coefficients_[k], q);
+    SymmetricMatrix<RT,embeddedDim> hesse = TargetSpace::secondDerivativeOfDistanceSquaredWRTSecondArgument(coefficients_[k], q);
     for (int i=0; i<embeddedDim; i++)
       for (int j=0; j<=i; j++)
         dqdwF(i, j, k) = dqdwF(j, i, k) = hesse(i,j);
@@ -838,5 +842,6 @@ private:
   std::unique_ptr<LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,Rotation<field_type,3> > > orientationFEFunction_;
 };
 
+}  // namespace Dune::GFE
 
 #endif
diff --git a/dune/gfe/functions/localprojectedfefunction.hh b/dune/gfe/functions/localprojectedfefunction.hh
index 936a8a78..c926e70c 100644
--- a/dune/gfe/functions/localprojectedfefunction.hh
+++ b/dune/gfe/functions/localprojectedfefunction.hh
@@ -13,9 +13,8 @@
 #include <dune/gfe/spaces/productmanifold.hh>
 #include <dune/gfe/spaces/rotation.hh>
 
-namespace Dune {
-
-  namespace GFE {
+namespace Dune::GFE
+{
 
     /** \brief Interpolate in an embedding Euclidean space, and project back onto the Riemannian manifold
      *
@@ -675,7 +674,6 @@ namespace Dune {
 
     };
 
-  }
+}  // namespace Dune::GFE
 
-}
 #endif
diff --git a/dune/gfe/functions/localquickanddirtyfefunction.hh b/dune/gfe/functions/localquickanddirtyfefunction.hh
index 91c3c39d..51d8003b 100644
--- a/dune/gfe/functions/localquickanddirtyfefunction.hh
+++ b/dune/gfe/functions/localquickanddirtyfefunction.hh
@@ -10,9 +10,8 @@
 #include <dune/gfe/linearalgebra.hh>
 #include <dune/gfe/spaces/rotation.hh>
 
-namespace Dune {
-
-  namespace GFE {
+namespace Dune::GFE
+{
 
     /** \brief Interpolate on a manifold, as fast as we can
      *
@@ -149,7 +148,6 @@ namespace Dune {
         return TargetSpace(c);
       }
     }
-  }   // namespace GFE
 
-}   // namespace Dune
+}   // namespace Dune::GFE
 #endif
diff --git a/dune/gfe/geodesicdifference.hh b/dune/gfe/geodesicdifference.hh
index 38788d1e..2afe61c1 100644
--- a/dune/gfe/geodesicdifference.hh
+++ b/dune/gfe/geodesicdifference.hh
@@ -5,6 +5,10 @@
 
 #include <dune/istl/bvector.hh>
 
+
+namespace Dune::GFE
+{
+
 template <class TargetSpace>
 Dune::BlockVector<typename TargetSpace::TangentVector> computeGeodesicDifference(const std::vector<TargetSpace>& a,
                                                                                  const std::vector<TargetSpace>& b)
@@ -24,4 +28,6 @@ Dune::BlockVector<typename TargetSpace::TangentVector> computeGeodesicDifference
   return result;
 }
 
+}  // namespace Dune::GFE
+
 #endif
diff --git a/dune/gfe/geodesicfefunctionadaptor.hh b/dune/gfe/geodesicfefunctionadaptor.hh
index e070f84d..e08c16a6 100644
--- a/dune/gfe/geodesicfefunctionadaptor.hh
+++ b/dune/gfe/geodesicfefunctionadaptor.hh
@@ -6,6 +6,10 @@
 
 #include "localgeodesicfefunction.hh"
 
+
+namespace Dune::GFE
+{
+
 template <class Basis, class TargetSpace>
 class GeodesicFEFunctionAdaptor
 {
@@ -208,4 +212,6 @@ public:
 
 };
 
+}  // namespace Dune::GFE
+
 #endif
diff --git a/dune/gfe/geometries/moebiusstrip.hh b/dune/gfe/geometries/moebiusstrip.hh
index 5d345c9d..509155b1 100644
--- a/dune/gfe/geometries/moebiusstrip.hh
+++ b/dune/gfe/geometries/moebiusstrip.hh
@@ -7,7 +7,8 @@
 #include <dune/common/fvector.hh>
 #include <dune/curvedgrid/gridfunctions/analyticgridfunction.hh>
 
-namespace Dune {
+namespace Dune::GFE
+{
 
   /// \brief Functor representing a Möbius strip in 3D, where the base circle is the circle in the x-y-plane with radius_ around 0,0,0
   template <class T = double>
@@ -101,6 +102,6 @@ namespace Dune {
     return analyticGridFunction<Grid>(MoebiusStripProjection<T>{radius});
   }
 
-} // end namespace Dune
+}  // namespace Dune::GFE
 
 #endif // DUNE_GFE_MOEBIUSSTRIP_GRIDFUNCTION_HH
diff --git a/dune/gfe/geometries/twistedstrip.hh b/dune/gfe/geometries/twistedstrip.hh
index 9d502aeb..06f01d3d 100644
--- a/dune/gfe/geometries/twistedstrip.hh
+++ b/dune/gfe/geometries/twistedstrip.hh
@@ -7,7 +7,8 @@
 #include <dune/common/fvector.hh>
 #include <dune/curvedgrid/gridfunctions/analyticgridfunction.hh>
 
-namespace Dune {
+namespace Dune::GFE
+{
 
   /// \brief Functor representing a twisted strip in 3D, with length length_ and nTwists twists
   template <class T = double>
@@ -99,6 +100,6 @@ namespace Dune {
     return analyticGridFunction<Grid>(TwistedStripProjection<T>{length, nTwists});
   }
 
-} // end namespace Dune
+}  // namespace Dune::GFE
 
 #endif // DUNE_GFE_TWISTEDSTRIP_GRIDFUNCTION_HH
diff --git a/dune/gfe/gfedifferencenormsquared.hh b/dune/gfe/gfedifferencenormsquared.hh
index a6d15123..9ab50398 100644
--- a/dune/gfe/gfedifferencenormsquared.hh
+++ b/dune/gfe/gfedifferencenormsquared.hh
@@ -15,6 +15,10 @@
 
 #include <dune/gfe/globalgeodesicfefunction.hh>
 
+
+namespace Dune::GFE
+{
+
 template <class Basis, class TargetSpace>
 class GFEDifferenceNormSquared {
 
@@ -285,4 +289,6 @@ public:
   }
 };
 
+}  // namespace Dune::GFE
+
 #endif
diff --git a/dune/gfe/gramschmidtsolver.hh b/dune/gfe/gramschmidtsolver.hh
index 69104265..d99f55c7 100644
--- a/dune/gfe/gramschmidtsolver.hh
+++ b/dune/gfe/gramschmidtsolver.hh
@@ -5,6 +5,10 @@
 
 #include <dune/gfe/symmetricmatrix.hh>
 
+
+namespace Dune::GFE
+{
+
 /** \brief Direct solver for a dense symmetric linear system, using an orthonormal basis
  *
  * This solver computes an A-orthonormal basis, and uses that to compute the solution
@@ -24,7 +28,7 @@ class GramSchmidtSolver
    * \param matrix The matrix inducing the matrix norm
    * \param[in,out] v The vector to normalize
    */
-  static void normalize(const Dune::SymmetricMatrix<field_type,embeddedDim>& matrix,
+  static void normalize(const SymmetricMatrix<field_type,embeddedDim>& matrix,
                         Dune::FieldVector<field_type,embeddedDim>& v)
   {
     using std::sqrt;
@@ -36,7 +40,7 @@ class GramSchmidtSolver
    *
    * \param matrix The matrix the defines the scalar product
    */
-  static void project(const Dune::SymmetricMatrix<field_type,embeddedDim>& matrix,
+  static void project(const SymmetricMatrix<field_type,embeddedDim>& matrix,
                       const Dune::FieldVector<field_type,embeddedDim>& vi,
                       Dune::FieldVector<field_type,embeddedDim>& vj)
   {
@@ -59,7 +63,7 @@ public:
    * \param basis Any basis of the orthogonal complement of the kernel,
    *              used as the input for the Gram-Schmidt orthogonalization process
    */
-  GramSchmidtSolver(const Dune::SymmetricMatrix<field_type,embeddedDim>& matrix,
+  GramSchmidtSolver(const SymmetricMatrix<field_type,embeddedDim>& matrix,
                     const Dune::FieldMatrix<field_type,rank,embeddedDim>& basis)
     : orthonormalBasis_(basis)
   {
@@ -95,7 +99,7 @@ public:
 
    * \param basis Any basis of the space, used as the input for the Gram-Schmidt orthogonalization process
    */
-  static void solve(const Dune::SymmetricMatrix<field_type,embeddedDim>& matrix,
+  static void solve(const SymmetricMatrix<field_type,embeddedDim>& matrix,
                     Dune::FieldVector<field_type,embeddedDim>& x,
                     const Dune::FieldVector<field_type,embeddedDim>& rhs,
                     const Dune::FieldMatrix<field_type,rank,embeddedDim>& basis)
@@ -131,4 +135,6 @@ private:
 
 };
 
+}  // namespace Dune::GFE
+
 #endif
diff --git a/dune/gfe/linearalgebra.hh b/dune/gfe/linearalgebra.hh
index 20c1836b..486ea0a0 100644
--- a/dune/gfe/linearalgebra.hh
+++ b/dune/gfe/linearalgebra.hh
@@ -14,9 +14,8 @@
 //  Various matrix methods
 ///////////////////////////////////////////////////////////////////////////////////////////
 
-namespace Dune {
-
-  namespace GFE {
+namespace Dune::GFE
+{
 
 #if ADOLC_ADOUBLE_H
     /** \brief Calculates ret = s*A, where A has as field_type of adouble.
@@ -231,8 +230,7 @@ namespace Dune {
       std::generate(vec.begin(), vec.end(), rand);
       return vec;
     }
-  }
-}
 
+}  // namespace Dune::GFE
 
 #endif
diff --git a/dune/gfe/maxnormtrustregion.hh b/dune/gfe/maxnormtrustregion.hh
index 50d649da..a4f11d21 100644
--- a/dune/gfe/maxnormtrustregion.hh
+++ b/dune/gfe/maxnormtrustregion.hh
@@ -5,6 +5,10 @@
 
 #include <dune/solvers/common/boxconstraint.hh>
 
+
+namespace Dune::GFE
+{
+
 template <int blocksize, class field_type=double>
 class MaxNormTrustRegion
 {
@@ -96,4 +100,6 @@ private:
 
 };
 
+}  // namespace Dune::GFE
+
 #endif
diff --git a/dune/gfe/mixedriemanniantrsolver.cc b/dune/gfe/mixedriemanniantrsolver.cc
index e1170435..182a0640 100644
--- a/dune/gfe/mixedriemanniantrsolver.cc
+++ b/dune/gfe/mixedriemanniantrsolver.cc
@@ -25,7 +25,7 @@ template <class GridType,
     class Basis,
     class Basis0, class TargetSpace0,
     class Basis1, class TargetSpace1>
-void MixedRiemannianTrustRegionSolver<GridType,Basis,Basis0,TargetSpace0,Basis1,TargetSpace1>::
+void Dune::GFE::MixedRiemannianTrustRegionSolver<GridType,Basis,Basis0,TargetSpace0,Basis1,TargetSpace1>::
 setup(const GridType& grid,
       const MixedGFEAssembler<Basis, TargetSpace>* assembler,
       const Basis0& tmpBasis0,
@@ -262,7 +262,7 @@ template <class GridType,
     class Basis,
     class Basis0, class TargetSpace0,
     class Basis1, class TargetSpace1>
-void MixedRiemannianTrustRegionSolver<GridType,Basis,Basis0,TargetSpace0,Basis1,TargetSpace1>::solve()
+void Dune::GFE::MixedRiemannianTrustRegionSolver<GridType,Basis,Basis0,TargetSpace0,Basis1,TargetSpace1>::solve()
 {
   int argc = 0;
   char** argv;
diff --git a/dune/gfe/mixedriemanniantrsolver.hh b/dune/gfe/mixedriemanniantrsolver.hh
index 2c17bbcb..4a1ecd61 100644
--- a/dune/gfe/mixedriemanniantrsolver.hh
+++ b/dune/gfe/mixedriemanniantrsolver.hh
@@ -22,6 +22,10 @@
 #include <dune/gfe/assemblers/mixedgfeassembler.hh>
 #include <dune/gfe/spaces/productmanifold.hh>
 
+
+namespace Dune::GFE
+{
+
 /** \brief Riemannian trust-region solver for geodesic finite-element problems */
 template <class GridType,
     class Basis,
@@ -185,6 +189,8 @@ protected:
   Statistics statistics_;
 };
 
+}  // namespace Dune::GFE
+
 #include "mixedriemanniantrsolver.cc"
 
 #endif
diff --git a/dune/gfe/neumannenergy.hh b/dune/gfe/neumannenergy.hh
index de6c78a6..c10fa104 100644
--- a/dune/gfe/neumannenergy.hh
+++ b/dune/gfe/neumannenergy.hh
@@ -8,7 +8,9 @@
 #include <dune/gfe/assemblers/localenergy.hh>
 #include <dune/gfe/spaces/productmanifold.hh>
 
-namespace Dune::GFE {
+namespace Dune::GFE
+{
+
   /** \brief Integrate a density over a part of the domain boundary
    *
    * This is typically used to implement Neumann boundary conditions for Cosserat materials.
diff --git a/dune/gfe/parallel/globalmapper.hh b/dune/gfe/parallel/globalmapper.hh
index 9aad9c30..503728af 100644
--- a/dune/gfe/parallel/globalmapper.hh
+++ b/dune/gfe/parallel/globalmapper.hh
@@ -7,7 +7,8 @@
 #include <dune/parmg/parallel/datahandle.hh>
 #endif
 
-namespace Dune {
+namespace Dune::GFE
+{
 
   template <class Basis>
   class GlobalMapper
@@ -63,5 +64,6 @@ namespace Dune {
     std::size_t size_;
   };
 
-}
+}  // namespace Dune::GFE
+
 #endif   // DUNE_GFE_PARALLEL_GLOBALMAPPER_HH
diff --git a/dune/gfe/parallel/globalp1mapper.hh b/dune/gfe/parallel/globalp1mapper.hh
index 90efa155..7b83f634 100644
--- a/dune/gfe/parallel/globalp1mapper.hh
+++ b/dune/gfe/parallel/globalp1mapper.hh
@@ -6,7 +6,8 @@
 
 #include <dune/gfe/parallel/globalmapper.hh>
 
-namespace Dune {
+namespace Dune::GFE
+{
 
   template <class Basis>
   class GlobalP1Mapper : public GlobalMapper<Basis>
@@ -79,5 +80,6 @@ namespace Dune {
 #endif
   };
 
-}
+}  // namespace Dune::GFE
+
 #endif /* DUNE_GFE_PARALLEL_GLOBALP1MAPPER_HH */
diff --git a/dune/gfe/parallel/globalp2mapper.hh b/dune/gfe/parallel/globalp2mapper.hh
index 3bc1fd99..7dfba294 100644
--- a/dune/gfe/parallel/globalp2mapper.hh
+++ b/dune/gfe/parallel/globalp2mapper.hh
@@ -7,10 +7,11 @@
 
 #include <dune/gfe/parallel/globalmapper.hh>
 
-namespace Dune {
+namespace Dune::GFE
+{
 
   template <class Basis>
-  class GlobalP2Mapper : public Dune::GlobalMapper<Basis>
+  class GlobalP2Mapper : public GlobalMapper<Basis>
   {
     using GridView = typename Basis::GridView;
     using P2BasisMapper = Functions::LagrangeBasis<GridView,2>;
@@ -119,5 +120,7 @@ namespace Dune {
 
     P2BasisMapper p2Mapper_;
   };
-}
+
+}  // namespace Dune::GFE
+
 #endif   // DUNE_GFE_PARALLEL_GLOBALP2MAPPER_HH
diff --git a/dune/gfe/parallel/mapperfactory.hh b/dune/gfe/parallel/mapperfactory.hh
index 4bd2bf23..345f527a 100644
--- a/dune/gfe/parallel/mapperfactory.hh
+++ b/dune/gfe/parallel/mapperfactory.hh
@@ -19,7 +19,7 @@ namespace Dune::GFE
   template <typename GridView>
   struct MapperFactory<Functions::LagrangeBasis<GridView,1> >
   {
-    typedef Dune::GlobalP1Mapper<Functions::LagrangeBasis<GridView,1> > GlobalMapper;
+    typedef GlobalP1Mapper<Functions::LagrangeBasis<GridView,1> > GlobalMapper;
     typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView> LocalMapper;
     static LocalMapper createLocalMapper(const GridView& gridView)
     {
@@ -30,7 +30,7 @@ namespace Dune::GFE
   template <typename GridView>
   struct MapperFactory<Dune::Functions::LagrangeBasis<GridView,2> >
   {
-    typedef Dune::GlobalP2Mapper<Functions::LagrangeBasis<GridView,2> > GlobalMapper;
+    typedef GlobalP2Mapper<Functions::LagrangeBasis<GridView,2> > GlobalMapper;
     typedef P2BasisMapper<GridView> LocalMapper;
     static LocalMapper createLocalMapper(const GridView& gridView)
     {
diff --git a/dune/gfe/parallel/matrixcommunicator.hh b/dune/gfe/parallel/matrixcommunicator.hh
index 920f066f..18289e1f 100644
--- a/dune/gfe/parallel/matrixcommunicator.hh
+++ b/dune/gfe/parallel/matrixcommunicator.hh
@@ -9,6 +9,9 @@
 #include <dune/gfe/parallel/mpifunctions.hh>
 
 
+namespace Dune::GFE
+{
+
 template<typename RowGlobalMapper, typename GridView1, typename GridView2, typename MatrixType, typename LocalMapper1, typename LocalMapper2, typename ColumnGlobalMapper=RowGlobalMapper>
 class MatrixCommunicator {
 
@@ -174,4 +177,6 @@ private:
   std::vector<TransferMatrixTuple> globalMatrixEntries;
 };
 
+}  // namespace Dune::GFE
+
 #endif
diff --git a/dune/gfe/parallel/mpifunctions.hh b/dune/gfe/parallel/mpifunctions.hh
index 3092e036..c5854d8c 100644
--- a/dune/gfe/parallel/mpifunctions.hh
+++ b/dune/gfe/parallel/mpifunctions.hh
@@ -4,6 +4,10 @@
 #include <numeric>
 #include <vector>
 
+
+namespace Dune::GFE
+{
+
 struct MPIFunctions {
   static std::vector<int> offsetsFromSizes(const std::vector<int>& sizes) {
     std::vector<int> offsets(sizes.size());
@@ -66,5 +70,6 @@ struct MPIFunctions {
   }
 };
 
+}  // namespace Dune::GFE
 
 #endif
diff --git a/dune/gfe/parallel/p2mapper.hh b/dune/gfe/parallel/p2mapper.hh
index 5dcf6969..34098d7f 100644
--- a/dune/gfe/parallel/p2mapper.hh
+++ b/dune/gfe/parallel/p2mapper.hh
@@ -8,6 +8,10 @@
 
 #include <dune/functions/functionspacebases/lagrangebasis.hh>
 
+
+namespace Dune::GFE
+{
+
 /** \brief Mimic a dune-grid mapper for a P2 space, using the dune-functions dof ordering of such a space
  */
 template<class GridView>
@@ -49,4 +53,6 @@ public:
   Dune::Functions::LagrangeBasis<GridView,2> p2Basis_;
 };
 
+}  // namespace Dune::GFE
+
 #endif
diff --git a/dune/gfe/parallel/vectorcommunicator.hh b/dune/gfe/parallel/vectorcommunicator.hh
index d83cc649..b38ed711 100644
--- a/dune/gfe/parallel/vectorcommunicator.hh
+++ b/dune/gfe/parallel/vectorcommunicator.hh
@@ -7,6 +7,9 @@
 #include <dune/gfe/parallel/mpifunctions.hh>
 
 
+namespace Dune::GFE
+{
+
 template<typename GUIndex, typename Communicator, typename VectorType>
 class VectorCommunicator {
 
@@ -101,4 +104,6 @@ private:
   std::vector<TransferVectorTuple> globalVectorEntries;
 };
 
+}  // namespace Dune::GFE
+
 #endif
diff --git a/dune/gfe/polardecomposition.hh b/dune/gfe/polardecomposition.hh
index b25e5416..6e186032 100644
--- a/dune/gfe/polardecomposition.hh
+++ b/dune/gfe/polardecomposition.hh
@@ -13,7 +13,8 @@
 //  Two Methods to compute the Polar Factor of a 3x3 matrix
 ///////////////////////////////////////////////////////////////////////////////////////////
 
-namespace Dune::GFE {
+namespace Dune::GFE
+{
 
   class PolarDecomposition
   {
@@ -321,6 +322,7 @@ namespace Dune::GFE {
       return v;
     }
   };
-}
+
+}  // namespace Dune::GFE
 
 #endif
diff --git a/dune/gfe/riemannianpnsolver.cc b/dune/gfe/riemannianpnsolver.cc
index 5bef95d7..2cee6d0e 100644
--- a/dune/gfe/riemannianpnsolver.cc
+++ b/dune/gfe/riemannianpnsolver.cc
@@ -19,7 +19,7 @@
 #include <dune/gfe/parallel/vectorcommunicator.hh>
 
 template <class Basis, class TargetSpace, class Assembler>
-void RiemannianProximalNewtonSolver<Basis, TargetSpace, Assembler>::
+void Dune::GFE::RiemannianProximalNewtonSolver<Basis, TargetSpace, Assembler>::
 setup(const GridType& grid,
       const Assembler* assembler,
       const SolutionType& x,
@@ -65,7 +65,7 @@ setup(const GridType& grid,
 }
 
 template <class Basis, class TargetSpace, class Assembler>
-void RiemannianProximalNewtonSolver<Basis,TargetSpace,Assembler>::
+void Dune::GFE::RiemannianProximalNewtonSolver<Basis,TargetSpace,Assembler>::
 setup(const GridType& grid,
       const Assembler* assembler,
       const SolutionType& x,
@@ -153,8 +153,8 @@ setup(const GridType& grid,
 
   // Write all intermediate solutions, if requested
   if (instrumented_
-      && dynamic_cast<IterativeSolver<CorrectionType>*>(innerSolver_.get()))
-    dynamic_cast<IterativeSolver<CorrectionType>*>(innerSolver_.get())->historyBuffer_ = instrumentedPath_ + "/mgHistory";
+      && dynamic_cast<::IterativeSolver<CorrectionType>*>(innerSolver_.get()))
+    dynamic_cast<::IterativeSolver<CorrectionType>*>(innerSolver_.get())->historyBuffer_ = instrumentedPath_ + "/mgHistory";
 
   // ////////////////////////////////////////////////////////////
   //    Create Hessian matrix and its occupation structure
@@ -168,7 +168,7 @@ setup(const GridType& grid,
 
 
 template <class Basis, class TargetSpace, class Assembler>
-void RiemannianProximalNewtonSolver<Basis,TargetSpace,Assembler>::solve()
+void Dune::GFE::RiemannianProximalNewtonSolver<Basis,TargetSpace,Assembler>::solve()
 {
   int rank = grid_->comm().rank();
 
diff --git a/dune/gfe/riemannianpnsolver.hh b/dune/gfe/riemannianpnsolver.hh
index 8cea9149..b1f4cd22 100644
--- a/dune/gfe/riemannianpnsolver.hh
+++ b/dune/gfe/riemannianpnsolver.hh
@@ -20,10 +20,13 @@
 #include <dune/gfe/parallel/mapperfactory.hh>
 
 
+namespace Dune::GFE
+{
+
 /** \brief Riemannian proximal-newton solver for geodesic finite-element problems */
 template <class Basis, class TargetSpace, class Assembler = GeodesicFEAssembler<Basis,TargetSpace> >
 class RiemannianProximalNewtonSolver
-  : public IterativeSolver<std::vector<TargetSpace>,
+  : public ::IterativeSolver<std::vector<TargetSpace>,
         Dune::BitSetVector<TargetSpace::TangentVector::dimension> >
 {
   typedef typename Basis::GridView::Grid GridType;
@@ -59,7 +62,7 @@ class RiemannianProximalNewtonSolver
 public:
 
   RiemannianProximalNewtonSolver()
-    : IterativeSolver<std::vector<TargetSpace>, Dune::BitSetVector<blocksize> >(0,100,NumProc::FULL),
+    : ::IterativeSolver<std::vector<TargetSpace>, Dune::BitSetVector<blocksize> >(0,100,NumProc::FULL),
     hessianMatrix_(nullptr), h1SemiNorm_(NULL)
   {
     std::fill(scaling_.begin(), scaling_.end(), 1.0);
@@ -165,6 +168,8 @@ protected:
 
 };
 
+}  // namespace Dune::GFE
+
 #include "riemannianpnsolver.cc"
 
 #endif
diff --git a/dune/gfe/riemanniantrsolver.cc b/dune/gfe/riemanniantrsolver.cc
index 2b6e0c5b..92f5b9ed 100644
--- a/dune/gfe/riemanniantrsolver.cc
+++ b/dune/gfe/riemanniantrsolver.cc
@@ -27,7 +27,7 @@
 #include <dune/gfe/parallel/vectorcommunicator.hh>
 
 template <class Basis, class TargetSpace, class Assembler>
-void RiemannianTrustRegionSolver<Basis, TargetSpace, Assembler>::
+void Dune::GFE::RiemannianTrustRegionSolver<Basis, TargetSpace, Assembler>::
 setup(const GridType& grid,
       const Assembler* assembler,
       const SolutionType& x,
@@ -66,7 +66,7 @@ setup(const GridType& grid,
 }
 
 template <class Basis, class TargetSpace, class Assembler>
-void RiemannianTrustRegionSolver<Basis,TargetSpace,Assembler>::
+void Dune::GFE::RiemannianTrustRegionSolver<Basis,TargetSpace,Assembler>::
 setup(const GridType& grid,
       const Assembler* assembler,
       const SolutionType& x,
@@ -121,7 +121,7 @@ setup(const GridType& grid,
   // Hack: the two-norm may not scale all that well, but it is fast!
   auto baseNorm = std::make_shared<TwoNorm<CorrectionType> >();
 
-  auto baseSolver = std::make_shared<::LoopSolver<CorrectionType> >(baseSolverStep,
+  auto baseSolver = std::make_shared<Solvers::LoopSolver<CorrectionType> >(baseSolverStep,
                                                                     baseIterations,
                                                                     baseTolerance,
                                                                     baseNorm,
@@ -179,7 +179,7 @@ setup(const GridType& grid,
 #endif
   h1SemiNorm_ = std::make_shared<H1SemiNorm<CorrectionType> >(A);
 
-  innerSolver_ = std::make_shared<::LoopSolver<CorrectionType> >(mmgStep,
+  innerSolver_ = std::make_shared<Solvers::LoopSolver<CorrectionType> >(mmgStep,
                                                                  innerIterations_,
                                                                  innerTolerance_,
                                                                  h1SemiNorm_,
@@ -204,8 +204,8 @@ setup(const GridType& grid,
 
   // Write all intermediate solutions, if requested
   if (instrumented_
-      && dynamic_cast<IterativeSolver<CorrectionType>*>(innerSolver_.get()))
-    dynamic_cast<IterativeSolver<CorrectionType>*>(innerSolver_.get())->historyBuffer_ = instrumentedPath_ + "/mgHistory";
+      && dynamic_cast<Solvers::IterativeSolver<CorrectionType>*>(innerSolver_.get()))
+    dynamic_cast<Solvers::IterativeSolver<CorrectionType>*>(innerSolver_.get())->historyBuffer_ = instrumentedPath_ + "/mgHistory";
 
   // ////////////////////////////////////////////////////////////
   //    Create Hessian matrix and its occupation structure
@@ -248,7 +248,7 @@ setup(const GridType& grid,
 #if HAVE_MPI
     // If we are on more than 1 processors, join all local transfer matrices on rank 0,
     // and construct a single global transfer operator there.
-    typedef Dune::GlobalP1Mapper<Dune::Functions::LagrangeBasis<typename Basis::GridView,1> > GlobalLeafP1Mapper;
+    typedef GFE::GlobalP1Mapper<Dune::Functions::LagrangeBasis<typename Basis::GridView,1> > GlobalLeafP1Mapper;
     GlobalLeafP1Mapper p1Index(grid_->leafGridView());
 
     typedef Dune::MultipleCodimMultipleGeomTypeMapper<typename GridType::LeafGridView> LeafP1LocalMapper;
@@ -282,7 +282,7 @@ setup(const GridType& grid,
     // If we are on more than 1 processors, join all local transfer matrices on rank 0,
     // and construct a single global transfer operator there.
     typedef Dune::Functions::LagrangeBasis<typename GridType::LevelGridView, 1> FEBasis;
-    typedef Dune::GlobalP1Mapper<FEBasis> GlobalLevelP1Mapper;
+    typedef GFE::GlobalP1Mapper<FEBasis> GlobalLevelP1Mapper;
     GlobalLevelP1Mapper fineGUIndex(grid_->levelGridView(i+1));
     GlobalLevelP1Mapper coarseGUIndex(grid_->levelGridView(i));
 
@@ -328,15 +328,15 @@ setup(const GridType& grid,
 
 
 template <class Basis, class TargetSpace, class Assembler>
-void RiemannianTrustRegionSolver<Basis,TargetSpace,Assembler>::solve()
+void Dune::GFE::RiemannianTrustRegionSolver<Basis,TargetSpace,Assembler>::solve()
 {
   int rank = grid_->comm().rank();
 
   MonotoneMGStep<MatrixType,CorrectionType>* mgStep = nullptr;    // Non-shared pointer -- the innerSolver keeps the ownership
 
   // if the inner solver is a monotone multigrid set up a max-norm trust-region
-  if (dynamic_cast<LoopSolver<CorrectionType>*>(innerSolver_.get())) {
-    auto loopSolver = std::dynamic_pointer_cast<LoopSolver<CorrectionType> >(innerSolver_);
+  if (dynamic_cast<Solvers::LoopSolver<CorrectionType>*>(innerSolver_.get())) {
+    auto loopSolver = std::dynamic_pointer_cast<Solvers::LoopSolver<CorrectionType> >(innerSolver_);
     mgStep = dynamic_cast<MonotoneMGStep<MatrixType,CorrectionType>*>(&loopSolver->getIterationStep());
   }
 
diff --git a/dune/gfe/riemanniantrsolver.hh b/dune/gfe/riemanniantrsolver.hh
index d2b000ac..7c122f3e 100644
--- a/dune/gfe/riemanniantrsolver.hh
+++ b/dune/gfe/riemanniantrsolver.hh
@@ -20,10 +20,13 @@
 #include <dune/gfe/parallel/mapperfactory.hh>
 
 
+namespace Dune::GFE
+{
+
 /** \brief Riemannian trust-region solver for geodesic finite-element problems */
 template <class Basis, class TargetSpace, class Assembler = GeodesicFEAssembler<Basis,TargetSpace> >
 class RiemannianTrustRegionSolver
-  : public IterativeSolver<std::vector<TargetSpace>,
+  : public ::IterativeSolver<std::vector<TargetSpace>,
         Dune::BitSetVector<TargetSpace::TangentVector::dimension> >
 {
   typedef typename Basis::GridView::Grid GridType;
@@ -59,7 +62,7 @@ class RiemannianTrustRegionSolver
 public:
 
   RiemannianTrustRegionSolver()
-    : IterativeSolver<std::vector<TargetSpace>, Dune::BitSetVector<blocksize> >(0,100,NumProc::FULL),
+    : ::IterativeSolver<std::vector<TargetSpace>, Dune::BitSetVector<blocksize> >(0,100,NumProc::FULL),
     hessianMatrix_(nullptr), h1SemiNorm_(NULL)
   {
     std::fill(scaling_.begin(), scaling_.end(), 1.0);
@@ -179,6 +182,8 @@ protected:
 
 };
 
+}  // namespace Dune::GFE
+
 #include "riemanniantrsolver.cc"
 
 #endif
diff --git a/dune/gfe/rodfactory.hh b/dune/gfe/rodfactory.hh
index bea34db0..fc7c516d 100644
--- a/dune/gfe/rodfactory.hh
+++ b/dune/gfe/rodfactory.hh
@@ -15,6 +15,10 @@
 #include <dune/gfe/spaces/realtuple.hh>
 #include <dune/gfe/spaces/rotation.hh>
 
+
+namespace Dune::GFE
+{
+
 /** \brief A factory class that implements various ways to create rod configurations
  */
 
@@ -233,4 +237,6 @@ private:
   const GridView gridView_;
 };
 
+}  // namespace Dune::GFE
+
 #endif
diff --git a/dune/gfe/skewmatrix.hh b/dune/gfe/skewmatrix.hh
index c45484b5..6f3305f2 100644
--- a/dune/gfe/skewmatrix.hh
+++ b/dune/gfe/skewmatrix.hh
@@ -1,6 +1,9 @@
 #ifndef DUNE_SKEW_MATRIX_HH
 #define DUNE_SKEW_MATRIX_HH
 
+namespace Dune::GFE
+{
+
 /** \brief Static dense skew-symmetric matrix */
 template <class T, int N>
 class SkewMatrix
@@ -98,4 +101,6 @@ private:
 
 };
 
+}  // namespace Dune::GFE
+
 #endif
diff --git a/dune/gfe/spaces/hyperbolichalfspacepoint.hh b/dune/gfe/spaces/hyperbolichalfspacepoint.hh
index e4b21e6a..ab841d09 100644
--- a/dune/gfe/spaces/hyperbolichalfspacepoint.hh
+++ b/dune/gfe/spaces/hyperbolichalfspacepoint.hh
@@ -9,6 +9,10 @@
 
 #include <dune/gfe/tensor3.hh>
 
+
+namespace Dune::GFE
+{
+
 /** \brief A point in the hyperbolic half-space H^N
 
     \tparam N Dimension of the hyperbolic half-space
@@ -549,4 +553,6 @@ private:
   Dune::FieldVector<T,N> data_;
 };
 
+}  // namespace Dune::GFE
+
 #endif
diff --git a/dune/gfe/spaces/orthogonalmatrix.hh b/dune/gfe/spaces/orthogonalmatrix.hh
index 71a8f76e..145992f7 100644
--- a/dune/gfe/spaces/orthogonalmatrix.hh
+++ b/dune/gfe/spaces/orthogonalmatrix.hh
@@ -5,6 +5,10 @@
 
 #include <dune/gfe/skewmatrix.hh>
 
+
+namespace Dune::GFE
+{
+
 /** \brief An orthogonal \f$ n \times n \f$ matrix
  * \tparam T Type of the matrix entries
  * \tparam N Space dimension
@@ -177,4 +181,6 @@ private:
 
 };
 
+}  // namespace Dune::GFE
+
 #endif  // DUNE_GFE_SPACES_ORTHOGONALMATRIX_HH
diff --git a/dune/gfe/spaces/productmanifold.hh b/dune/gfe/spaces/productmanifold.hh
index 5bb4ef91..2f403295 100644
--- a/dune/gfe/spaces/productmanifold.hh
+++ b/dune/gfe/spaces/productmanifold.hh
@@ -208,7 +208,7 @@ namespace Dune::GFE
     static auto secondDerivativeOfDistanceSquaredWRTSecondArgument(const ProductManifold& a,
                                                                    const ProductManifold& b)
     {
-      Dune::SymmetricMatrix<field_type,embeddedDim> result;
+      SymmetricMatrix<field_type,embeddedDim> result;
       auto secDerivOfDistSqWRTSecArgFunctor = [] (auto& argsTuple, std::array<std::size_t,2>& posHelper, const auto& manifoldInt)
                                               {
                                                 auto& deriv   = std::get<0>(argsTuple);
diff --git a/dune/gfe/spaces/quaternion.hh b/dune/gfe/spaces/quaternion.hh
index 141ac08d..5aa8d2d6 100644
--- a/dune/gfe/spaces/quaternion.hh
+++ b/dune/gfe/spaces/quaternion.hh
@@ -7,6 +7,10 @@
 
 #include <dune/gfe/tensor3.hh>
 
+
+namespace Dune::GFE
+{
+
 template <class T>
 class Quaternion : public Dune::FieldVector<T,4>
 {
@@ -161,11 +165,13 @@ public:
 
 };
 
+}  // namespace Dune::GFE
+
 namespace Dune
 {
   /** \brief Specizalization needed to allow certain forms of matrix--quaternion multiplications */
   template< class T >
-  struct FieldTraits< Quaternion<T> >
+  struct FieldTraits< Dune::GFE::Quaternion<T> >
   {
     typedef typename FieldTraits<T>::field_type field_type;
     typedef typename FieldTraits<T>::real_type real_type;
diff --git a/dune/gfe/spaces/realtuple.hh b/dune/gfe/spaces/realtuple.hh
index 92a674e3..d66908ab 100644
--- a/dune/gfe/spaces/realtuple.hh
+++ b/dune/gfe/spaces/realtuple.hh
@@ -9,6 +9,9 @@
 #include <dune/gfe/symmetricmatrix.hh>
 
 
+namespace Dune::GFE
+{
+
 /** \brief Implement a tuple of real numbers as a Riemannian manifold
 
    Currently this class only exists for testing purposes.
@@ -145,9 +148,9 @@ public:
 
      Unlike the distance itself the squared distance is differentiable at zero
    */
-  static Dune::SymmetricMatrix<T,N> secondDerivativeOfDistanceSquaredWRTSecondArgument(const RealTuple& a, const RealTuple& b) {
+  static SymmetricMatrix<T,N> secondDerivativeOfDistanceSquaredWRTSecondArgument(const RealTuple& a, const RealTuple& b) {
 
-    Dune::SymmetricMatrix<T,N> result;
+    SymmetricMatrix<T,N> result;
     for (int i=0; i<N; i++)
       for (int j=0; j<=i; j++)
         result(i,j) = 2*(i==j);
@@ -254,4 +257,6 @@ private:
 
 };
 
+}  // namespace Dune::GFE
+
 #endif
diff --git a/dune/gfe/spaces/rotation.hh b/dune/gfe/spaces/rotation.hh
index d5069d70..6a1f9d63 100644
--- a/dune/gfe/spaces/rotation.hh
+++ b/dune/gfe/spaces/rotation.hh
@@ -18,6 +18,10 @@
 #include <dune/gfe/spaces/quaternion.hh>
 #include <dune/gfe/spaces/unitvector.hh>
 
+
+namespace Dune::GFE
+{
+
 template <class T, int dim>
 class Rotation
 {};
@@ -772,13 +776,13 @@ public:
   }
 
   /** \brief Compute the Hessian of the squared distance function keeping the first argument fixed */
-  static Dune::SymmetricMatrix<T,4> secondDerivativeOfDistanceSquaredWRTSecondArgument(const Rotation<T,3>& p, const Rotation<T,3>& q) {
+  static SymmetricMatrix<T,4> secondDerivativeOfDistanceSquaredWRTSecondArgument(const Rotation<T,3>& p, const Rotation<T,3>& q) {
 
     T sp = p.globalCoordinates() * q.globalCoordinates();
 
     EmbeddedTangentVector pProjected = q.projectOntoTangentSpace(p.globalCoordinates());
 
-    Dune::SymmetricMatrix<T,4> A;
+    SymmetricMatrix<T,4> A;
     for (int i=0; i<4; i++)
       for (int j=0; j<=i; j++)
         A(i,j) = pProjected[i]*pProjected[j];
@@ -787,7 +791,7 @@ public:
     A *= 4*UnitVector<T,4>::secondDerivativeOfArcCosSquared(abs(sp));
 
     // Compute matrix B (see notes)
-    Dune::SymmetricMatrix<T,4> Pq;
+    SymmetricMatrix<T,4> Pq;
     for (int i=0; i<4; i++)
       for (int j=0; j<=i; j++)
         Pq(i,j) = (i==j) - q.globalCoordinates()[i]*q.globalCoordinates()[j];
@@ -1271,6 +1275,6 @@ public:
   };
 };
 
-
+}  // namespace Dune::GFE
 
 #endif
diff --git a/dune/gfe/spaces/unitvector.hh b/dune/gfe/spaces/unitvector.hh
index 31bd401e..68b7ed63 100644
--- a/dune/gfe/spaces/unitvector.hh
+++ b/dune/gfe/spaces/unitvector.hh
@@ -9,6 +9,10 @@
 #include <dune/gfe/tensor3.hh>
 #include <dune/gfe/symmetricmatrix.hh>
 
+
+namespace Dune::GFE
+{
+
 template <class T, int N>
 class Rotation;
 
@@ -291,13 +295,13 @@ public:
 
      Unlike the distance itself the squared distance is differentiable at zero
    */
-  static Dune::SymmetricMatrix<T,N> secondDerivativeOfDistanceSquaredWRTSecondArgument(const UnitVector& p, const UnitVector& q) {
+  static SymmetricMatrix<T,N> secondDerivativeOfDistanceSquaredWRTSecondArgument(const UnitVector& p, const UnitVector& q) {
 
     T sp = p.data_ * q.data_;
 
     Dune::FieldVector<T,N> pProjected = q.projectOntoTangentSpace(p.globalCoordinates());
 
-    Dune::SymmetricMatrix<T,N> A;
+    SymmetricMatrix<T,N> A;
     for (int i=0; i<N; i++)
       for (int j=0; j<=i; j++)
         A(i,j) = pProjected[i]*pProjected[j];
@@ -305,7 +309,7 @@ public:
     A *= secondDerivativeOfArcCosSquared(sp);
 
     // Compute matrix B (see notes)
-    Dune::SymmetricMatrix<T,N> Pq;
+    SymmetricMatrix<T,N> Pq;
     for (int i=0; i<N; i++)
       for (int j=0; j<=i; j++)
         Pq(i,j) = (i==j) - q.data_[i]*q.data_[j];
@@ -533,4 +537,6 @@ private:
   Dune::FieldVector<T,N> data_;
 };
 
+}  // namespace Dune::GFE
+
 #endif
diff --git a/dune/gfe/svd.hh b/dune/gfe/svd.hh
index 711093d2..14733672 100644
--- a/dune/gfe/svd.hh
+++ b/dune/gfe/svd.hh
@@ -12,6 +12,9 @@
 // int SIGN(const T& a, const T& b) {return 1;}
 #define SIGN(a,b)((b)>=0.0 ? fabs(a) : -fabs(a))
 
+namespace Dune::GFE
+{
+
 /** Computes (a^2 + b^2 )1/2 without destructive underflow or overflow. */
 template <class T>
 T pythag(T a, T b)
@@ -261,4 +264,6 @@ void svdcmp(Dune::FieldMatrix<T,m,n>& a_, Dune::FieldVector<T,n>& w, Dune::Field
       v_[i][j] = v[i+1][j+1];
 }
 
+}  // namespace Dune::GFE
+
 #endif
diff --git a/dune/gfe/symmetricmatrix.hh b/dune/gfe/symmetricmatrix.hh
index 7bc860ec..e4655863 100644
--- a/dune/gfe/symmetricmatrix.hh
+++ b/dune/gfe/symmetricmatrix.hh
@@ -4,7 +4,8 @@
 #include <dune/common/fvector.hh>
 #include <dune/common/fmatrix.hh>
 
-namespace Dune {
+namespace Dune::GFE
+{
 
   /** \brief  A class implementing a symmetric matrix with compile-time size
    *
@@ -93,5 +94,6 @@ namespace Dune {
     Dune::FieldVector<T,N*(N+1)/2> data_;
   };
 
-}
+}  // namespace Dune::GFE
+
 #endif
diff --git a/dune/gfe/targetspacertrsolver.cc b/dune/gfe/targetspacertrsolver.cc
index 9d2bc556..a5ca05c9 100644
--- a/dune/gfe/targetspacertrsolver.cc
+++ b/dune/gfe/targetspacertrsolver.cc
@@ -9,7 +9,7 @@
 #include <dune/gfe/gramschmidtsolver.hh>
 
 template <class TargetSpace>
-void TargetSpaceRiemannianTRSolver<TargetSpace>::
+void Dune::GFE::TargetSpaceRiemannianTRSolver<TargetSpace>::
 setup(const AverageDistanceAssembler<TargetSpace>* assembler,
       const TargetSpace& x,
       double tolerance,
@@ -25,7 +25,7 @@ setup(const AverageDistanceAssembler<TargetSpace>* assembler,
 
 
 template <class TargetSpace>
-void TargetSpaceRiemannianTRSolver<TargetSpace>::solve()
+void Dune::GFE::TargetSpaceRiemannianTRSolver<TargetSpace>::solve()
 {
   assert(minNumberOfIterations_ > 0);
 
diff --git a/dune/gfe/targetspacertrsolver.hh b/dune/gfe/targetspacertrsolver.hh
index e07714ab..fe435c1c 100644
--- a/dune/gfe/targetspacertrsolver.hh
+++ b/dune/gfe/targetspacertrsolver.hh
@@ -5,6 +5,10 @@
 
 #include <dune/gfe/symmetricmatrix.hh>
 
+
+namespace Dune::GFE
+{
+
 /** \brief Riemannian trust-region solver for geodesic finite-element problems
    \tparam TargetSpace The manifold that our functions take values in
  */
@@ -18,7 +22,7 @@ class TargetSpaceRiemannianTRSolver
   // Centralize the field type here
   typedef typename TargetSpace::ctype field_type;
 
-  typedef Dune::SymmetricMatrix<field_type, embeddedBlocksize> MatrixType;
+  typedef SymmetricMatrix<field_type, embeddedBlocksize> MatrixType;
   typedef Dune::FieldVector<field_type, embeddedBlocksize>     CorrectionType;
 
 public:
@@ -62,6 +66,8 @@ protected:
 
 };
 
+}  // namespace Dune::GFE
+
 #include "targetspacertrsolver.cc"
 
 #endif
diff --git a/dune/gfe/tensor3.hh b/dune/gfe/tensor3.hh
index 5bbd7225..01725b6d 100644
--- a/dune/gfe/tensor3.hh
+++ b/dune/gfe/tensor3.hh
@@ -8,6 +8,10 @@
 #include <array>
 #include <dune/common/fmatrix.hh>
 
+
+namespace Dune::GFE
+{
+
 /** \brief A third-rank tensor
  */
 template <class T, int N1, int N2, int N3>
@@ -193,6 +197,6 @@ inline std::ostream& operator<< (std::ostream& s, const Tensor3<T,N1,N2,N3>& ten
   return s;
 }
 
-
+}  // namespace Dune::GFE
 
 #endif
diff --git a/dune/gfe/tensorssd.hh b/dune/gfe/tensorssd.hh
index 0a9be67c..65a9baec 100644
--- a/dune/gfe/tensorssd.hh
+++ b/dune/gfe/tensorssd.hh
@@ -10,6 +10,10 @@
 #include <dune/common/fmatrix.hh>
 #include <dune/istl/matrix.hh>
 
+
+namespace Dune::GFE
+{
+
 /** \brief A third-rank tensor with two static (SS) and one dynamic (D) dimension
  *
  * \tparam T Type of the entries
@@ -127,5 +131,6 @@ inline std::ostream& operator<< (std::ostream& s, const TensorSSD<T,N1,N2>& tens
   return s;
 }
 
+}  // namespace Dune::GFE
 
 #endif
diff --git a/dune/gfe/trustregionmmgbasesolver.hh b/dune/gfe/trustregionmmgbasesolver.hh
index 9001e709..f14d00c8 100644
--- a/dune/gfe/trustregionmmgbasesolver.hh
+++ b/dune/gfe/trustregionmmgbasesolver.hh
@@ -9,6 +9,9 @@
 #include <dune/solvers/solvers/iterativesolver.hh>
 #include <dune/solvers/solvers/quadraticipopt.hh>
 
+namespace Dune::GFE
+{
+
 /** \brief Base solver for a monotone multigrid solver when used as the inner solver in a trust region method
  *
  * Monotone multigrid methods solve constrained problems even on the coarsest level.  Therefore, the choice
@@ -196,5 +199,6 @@ void TrustRegionMMGBaseSolver<MatrixType, VectorType>::solve()
 
 }
 
+}  // namespace Dune::GFE
 
 #endif
diff --git a/dune/gfe/vtkfile.hh b/dune/gfe/vtkfile.hh
index 3bfcd1f3..e3909183 100644
--- a/dune/gfe/vtkfile.hh
+++ b/dune/gfe/vtkfile.hh
@@ -17,9 +17,8 @@
 // For parallel infrastructure stuff:
 #include <dune/grid/io/file/vtk.hh>
 
-namespace Dune {
-
-  namespace GFE {
+namespace Dune::GFE
+{
 
     /** \brief A class representing a VTK file, but independent from the Dune grid interface
      *
@@ -367,8 +366,6 @@ namespace Dune {
 
     };
 
-  }
-
-}
+}  // namespace Dune::GFE
 
 #endif
diff --git a/src/compute-disc-error.cc b/src/compute-disc-error.cc
index edb32aac..910c0508 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,GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> > >::value)
+  if (std::is_same<TargetSpace,GFE::ProductManifold<GFE::RealTuple<double,3>,GFE::Rotation<double,3> > >::value)
   {
     double deformationL2ErrorSquared = 0;
     double orientationL2ErrorSquared = 0;
@@ -264,10 +264,10 @@ void measureDiscreteEOC(const GridView gridView,
         // Compute error of the orientation degrees of freedom
         // We need to transform from quaternion coordinates to matrix coordinates first.
         FieldMatrix<double,3,3> referenceValue, numericalValue;
-        Rotation<double,3> referenceRotation(std::array<double,4>{refValue[3], refValue[4], refValue[5], refValue[6]});
+        GFE::Rotation<double,3> referenceRotation(std::array<double,4>{refValue[3], refValue[4], refValue[5], refValue[6]});
         referenceRotation.matrix(referenceValue);
 
-        Rotation<double,3> numericalRotation(std::array<double,4>{numValue[3], numValue[4], numValue[5], numValue[6]});
+        GFE::Rotation<double,3> numericalRotation(std::array<double,4>{numValue[3], numValue[4], numValue[5], numValue[6]});
         numericalRotation.matrix(numericalValue);
 
         auto orientationDiff = referenceValue - numericalValue;
@@ -278,11 +278,11 @@ void measureDiscreteEOC(const GridView gridView,
         auto numericalDerQuat = localNumericalDerivative(localPos);
 
         // Transform to matrix coordinates
-        Tensor3<double,3,3,4> derivativeQuaternionToMatrixRef = Rotation<double,3>::derivativeOfQuaternionToMatrix(FieldVector<double,4>{refValue[3], refValue[4], refValue[5], refValue[6]});
-        Tensor3<double,3,3,4> derivativeQuaternionToMatrixNum = Rotation<double,3>::derivativeOfQuaternionToMatrix(FieldVector<double,4>{numValue[3], numValue[4], numValue[5], numValue[6]});
+        GFE::Tensor3<double,3,3,4> derivativeQuaternionToMatrixRef = GFE::Rotation<double,3>::derivativeOfQuaternionToMatrix(FieldVector<double,4>{refValue[3], refValue[4], refValue[5], refValue[6]});
+        GFE::Tensor3<double,3,3,4> derivativeQuaternionToMatrixNum = GFE::Rotation<double,3>::derivativeOfQuaternionToMatrix(FieldVector<double,4>{numValue[3], numValue[4], numValue[5], numValue[6]});
 
-        Tensor3<double,3,3,dim> refDerivative(0);
-        Tensor3<double,3,3,dim> numDerivative(0);
+        GFE::Tensor3<double,3,3,dim> refDerivative(0);
+        GFE::Tensor3<double,3,3,dim> numDerivative(0);
 
         for (int i=0; i<3; i++)
           for (int j=0; j<3; j++)
@@ -309,7 +309,7 @@ void measureDiscreteEOC(const GridView gridView,
               << "H^1 error orientation: " << std::sqrt(orientationH1ErrorSquared)
               << std::endl;
   }
-  else if constexpr (std::is_same<TargetSpace,Rotation<double,3> >::value)
+  else if constexpr (std::is_same<TargetSpace,GFE::Rotation<double,3> >::value)
   {
     double l2ErrorSquared = 0;
     double h1ErrorSquared = 0;
@@ -334,11 +334,11 @@ void measureDiscreteEOC(const GridView gridView,
 
         FieldMatrix<double,3,3> referenceValue, numericalValue;
         auto refValue = localReferenceSolution(qp.position());
-        Rotation<double,3> referenceRotation(refValue);
+        GFE::Rotation<double,3> referenceRotation(refValue);
         referenceRotation.matrix(referenceValue);
 
         auto numValue = localNumericalSolution(localPos);
-        Rotation<double,3> numericalRotation(numValue);
+        GFE::Rotation<double,3> numericalRotation(numValue);
         numericalRotation.matrix(numericalValue);
 
         auto diff = referenceValue - numericalValue;
@@ -349,11 +349,11 @@ void measureDiscreteEOC(const GridView gridView,
         auto numericalDerQuat = localNumericalDerivative(localPos);
 
         // Transform to matrix coordinates
-        Tensor3<double,3,3,4> derivativeQuaternionToMatrixRef = Rotation<double,3>::derivativeOfQuaternionToMatrix(refValue);
-        Tensor3<double,3,3,4> derivativeQuaternionToMatrixNum = Rotation<double,3>::derivativeOfQuaternionToMatrix(numValue);
+        GFE::Tensor3<double,3,3,4> derivativeQuaternionToMatrixRef = GFE::Rotation<double,3>::derivativeOfQuaternionToMatrix(refValue);
+        GFE::Tensor3<double,3,3,4> derivativeQuaternionToMatrixNum = GFE::Rotation<double,3>::derivativeOfQuaternionToMatrix(numValue);
 
-        Tensor3<double,3,3,dim> refDerivative(0);
-        Tensor3<double,3,3,dim> numDerivative(0);
+        GFE::Tensor3<double,3,3,dim> refDerivative(0);
+        GFE::Tensor3<double,3,3,dim> numDerivative(0);
 
         for (int i=0; i<3; i++)
           for (int j=0; j<3; j++)
@@ -473,7 +473,7 @@ void measureAnalyticalEOC(const GridView gridView,
   // TODO: We need to use a type-erasure wrapper here
   // Only used if // parameterSet["interpolationMethod"] == "geodesic"
   auto numericalSolutionGeodesic = GFE::EmbeddedGlobalGFEFunction<FEBasis,
-      LocalGeodesicFEFunction<dim, double, typename FEBasis::LocalView::Tree::FiniteElement, TargetSpace>,
+      GFE::LocalGeodesicFEFunction<dim, double, typename FEBasis::LocalView::Tree::FiniteElement, TargetSpace>,
       TargetSpace> (feBasis, x);
   auto localNumericalSolutionGeodesic = localFunction(numericalSolutionGeodesic);
 
@@ -493,7 +493,7 @@ void measureAnalyticalEOC(const GridView gridView,
 
   if (parameterSet["interpolationMethod"] == "geodesic")
     numericalSolution = std::make_unique<GFE::EmbeddedGlobalGFEFunction<FEBasis,
-        LocalGeodesicFEFunction<dim, double, typename FEBasis::LocalView::Tree::FiniteElement, TargetSpace>,
+        GFE::LocalGeodesicFEFunction<dim, double, typename FEBasis::LocalView::Tree::FiniteElement, TargetSpace>,
         TargetSpace> > (feBasis, x);
 
   if (parameterSet["interpolationMethod"] == "projected")
@@ -508,7 +508,7 @@ void measureAnalyticalEOC(const GridView gridView,
   // Compute errors in the L2 norm and the h1 seminorm.
   // SO(3)-valued maps need special treatment, because they are stored as quaternions,
   // but the errors need to be computed in matrix space.
-  if constexpr (std::is_same<TargetSpace,Rotation<double,3> >::value)
+  if constexpr (std::is_same<TargetSpace,GFE::Rotation<double,3> >::value)
   {
     constexpr int blocksize = TargetSpace::CoordinateType::dimension;
 
@@ -560,11 +560,11 @@ void measureAnalyticalEOC(const GridView gridView,
 #endif
 
         // Get error in matrix space
-        Rotation<double,3> numRotation(numValue);
+        GFE::Rotation<double,3> numRotation(numValue);
         FieldMatrix<double,3,3> numValueMatrix;
         numRotation.matrix(numValueMatrix);
 
-        Rotation<double,3> refRotation(refValue);
+        GFE::Rotation<double,3> refRotation(refValue);
         FieldMatrix<double,3,3> refValueMatrix;
         refRotation.matrix(refValueMatrix);
 
@@ -599,11 +599,11 @@ void measureAnalyticalEOC(const GridView gridView,
 #endif
 
         // Transform into matrix space
-        Tensor3<double,3,3,4> derivativeQuaternionToMatrixNum = Rotation<double,3>::derivativeOfQuaternionToMatrix(numValue);
-        Tensor3<double,3,3,4> derivativeQuaternionToMatrixRef = Rotation<double,3>::derivativeOfQuaternionToMatrix(refValue);
+        GFE::Tensor3<double,3,3,4> derivativeQuaternionToMatrixNum = GFE::Rotation<double,3>::derivativeOfQuaternionToMatrix(numValue);
+        GFE::Tensor3<double,3,3,4> derivativeQuaternionToMatrixRef = GFE::Rotation<double,3>::derivativeOfQuaternionToMatrix(refValue);
 
-        Tensor3<double,3,3,dim> numDerivative(0);
-        Tensor3<double,3,3,dim> refDerivative(0);
+        GFE::Tensor3<double,3,3,dim> numDerivative(0);
+        GFE::Tensor3<double,3,3,dim> refDerivative(0);
 
         for (int i=0; i<3; i++)
           for (int j=0; j<3; j++)
@@ -821,12 +821,12 @@ int main (int argc, char *argv[]) try
   case 1 :
     if (targetSpace=="RealTuple")
     {
-      measureEOC<GridType,RealTuple<double,1> >(grid,
+      measureEOC<GridType,GFE::RealTuple<double,1> >(grid,
                                                 referenceGrid,
                                                 parameterSet);
     } else if (targetSpace=="UnitVector")
     {
-      measureEOC<GridType,UnitVector<double,1> >(grid,
+      measureEOC<GridType,GFE::UnitVector<double,1> >(grid,
                                                  referenceGrid,
                                                  parameterSet);
     } else
@@ -836,23 +836,23 @@ int main (int argc, char *argv[]) try
   case 2 :
     if (targetSpace=="RealTuple")
     {
-      measureEOC<GridType,RealTuple<double,2> >(grid,
+      measureEOC<GridType,GFE::RealTuple<double,2> >(grid,
                                                 referenceGrid,
                                                 parameterSet);
     } else if (targetSpace=="UnitVector")
     {
-      measureEOC<GridType,UnitVector<double,2> >(grid,
+      measureEOC<GridType,GFE::UnitVector<double,2> >(grid,
                                                  referenceGrid,
                                                  parameterSet);
 #if 0
     } else if (targetSpace=="Rotation")
     {
-      measureEOC<GridType,Rotation<double,2> >(grid,
+      measureEOC<GridType,GFE::Rotation<double,2> >(grid,
                                                referenceGrid,
                                                parameterSet);
     } else if (targetSpace=="RigidBodyMotion")
     {
-      measureEOC<GridType,GFE::ProductManifold<RealTuple<double,2>,Rotation<double,2> > >(grid,
+      measureEOC<GridType,GFE::ProductManifold<GFE::RealTuple<double,2>,GFE::Rotation<double,2> > >(grid,
                                                                                           referenceGrid,
                                                                                           parameterSet);
 #endif
@@ -863,22 +863,22 @@ int main (int argc, char *argv[]) try
   case 3 :
     if (targetSpace=="RealTuple")
     {
-      measureEOC<GridType,RealTuple<double,3> >(grid,
+      measureEOC<GridType,GFE::RealTuple<double,3> >(grid,
                                                 referenceGrid,
                                                 parameterSet);
     } else if (targetSpace=="UnitVector")
     {
-      measureEOC<GridType,UnitVector<double,3> >(grid,
+      measureEOC<GridType,GFE::UnitVector<double,3> >(grid,
                                                  referenceGrid,
                                                  parameterSet);
     } else if (targetSpace=="Rotation")
     {
-      measureEOC<GridType,Rotation<double,3> >(grid,
+      measureEOC<GridType,GFE::Rotation<double,3> >(grid,
                                                referenceGrid,
                                                parameterSet);
     } else if (targetSpace=="RigidBodyMotion")
     {
-      measureEOC<GridType,GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> > >(grid,
+      measureEOC<GridType,GFE::ProductManifold<GFE::RealTuple<double,3>,GFE::Rotation<double,3> > >(grid,
                                                                                           referenceGrid,
                                                                                           parameterSet);
     } else
@@ -888,18 +888,18 @@ int main (int argc, char *argv[]) try
   case 4 :
     if (targetSpace=="RealTuple")
     {
-      measureEOC<GridType,RealTuple<double,4> >(grid,
+      measureEOC<GridType,GFE::RealTuple<double,4> >(grid,
                                                 referenceGrid,
                                                 parameterSet);
     } else if (targetSpace=="UnitVector")
     {
-      measureEOC<GridType,UnitVector<double,4> >(grid,
+      measureEOC<GridType,GFE::UnitVector<double,4> >(grid,
                                                  referenceGrid,
                                                  parameterSet);
 #if 0
     } else if (targetSpace=="Rotation")
     {
-      measureEOC<GridType,Rotation<double,4> >(grid,
+      measureEOC<GridType,GFE::Rotation<double,4> >(grid,
                                                referenceGrid,
                                                parameterSet);
     } else if (targetSpace=="RigidBodyMotion")
diff --git a/src/cosserat-continuum.cc b/src/cosserat-continuum.cc
index 3acdbd63..89176b23 100644
--- a/src/cosserat-continuum.cc
+++ b/src/cosserat-continuum.cc
@@ -107,7 +107,7 @@ static_assert(displacementOrder==rotationOrder, "displacement and rotation order
 #endif
 
 // Image space of the geodesic fe functions
-using TargetSpace = GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> >;
+using TargetSpace = GFE::ProductManifold<GFE::RealTuple<double,3>,GFE::Rotation<double,3> >;
 
 // Method to construct a Cosserat energy that matches the grid dimension.
 // This cannot be done inside the 'main' method, because 'constexpr if' only works
@@ -135,7 +135,7 @@ auto createCosseratEnergy(const ParameterTree& materialParameters,
     using Element = typename GridView::template Codim<0>::Entity;
     auto density = std::make_shared<GFE::CosseratShellDensity<Element, adouble> >(materialParameters);
 
-    return std::make_shared<NonplanarCosseratShellEnergy<Basis, 3, adouble, decltype(creator)> >(density, &creator);
+    return std::make_shared<GFE::NonplanarCosseratShellEnergy<Basis, 3, adouble, decltype(creator)> >(density, &creator);
   }
   else
   {
@@ -203,8 +203,8 @@ int main (int argc, char *argv[]) try
   const std::string resultPath          = parameterSet.get("resultPath", "");
 
   using namespace Dune::Indices;
-  using SolutionType = TupleVector<std::vector<RealTuple<double,3> >,
-      std::vector<Rotation<double,3> > >;
+  using SolutionType = TupleVector<std::vector<GFE::RealTuple<double,3> >,
+      std::vector<GFE::Rotation<double,3> > >;
 
   // ///////////////////////////////////////
   //    Create the grid
@@ -268,7 +268,7 @@ int main (int argc, char *argv[]) try
 
   using namespace Dune::Functions::BasisFactory;
 
-  const int dimRotation = Rotation<double,3>::embeddedDim;
+  const int dimRotation = GFE::Rotation<double,3>::embeddedDim;
   auto compositeBasis = makeBasis(
     gridView,
     composite(
@@ -438,7 +438,7 @@ int main (int argc, char *argv[]) try
 #ifdef PROJECTED_INTERPOLATION
     using LocalInterpolationRule  = GFE::LocalProjectedFEFunction<dim, double, DeformationFEBasis::LocalView::Tree::FiniteElement, TargetSpace>;
 #else
-    using LocalInterpolationRule  = LocalGeodesicFEFunction<dim, double, DeformationFEBasis::LocalView::Tree::FiniteElement, TargetSpace>;
+    using LocalInterpolationRule  = GFE::LocalGeodesicFEFunction<dim, double, DeformationFEBasis::LocalView::Tree::FiniteElement, TargetSpace>;
 #endif
     GFE::EmbeddedGlobalGFEFunction<InitialBasis,LocalInterpolationRule,TargetSpace> initialFunction(initialBasis,initialIterate);
     auto powerBasis = makeBasis(
@@ -477,25 +477,25 @@ int main (int argc, char *argv[]) try
   auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,3> >(deformationPowerBasis, displacement);
 
 #ifdef PROJECTED_INTERPOLATION
-  using RotationInterpolationRule = GFE::LocalProjectedFEFunction<dim, double, OrientationFEBasis::LocalView::Tree::FiniteElement, Rotation<double,3> >;
+  using RotationInterpolationRule = GFE::LocalProjectedFEFunction<dim, double, OrientationFEBasis::LocalView::Tree::FiniteElement, GFE::Rotation<double,3> >;
 #else
-  using RotationInterpolationRule = LocalGeodesicFEFunction<dim, double, OrientationFEBasis::LocalView::Tree::FiniteElement, Rotation<double,3> >;
+  using RotationInterpolationRule = GFE::LocalGeodesicFEFunction<dim, double, OrientationFEBasis::LocalView::Tree::FiniteElement, GFE::Rotation<double,3> >;
 #endif
 
-  GFE::EmbeddedGlobalGFEFunction<OrientationFEBasis, RotationInterpolationRule,Rotation<double,3> > orientationFunction(orientationFEBasis,
+  GFE::EmbeddedGlobalGFEFunction<OrientationFEBasis, RotationInterpolationRule,GFE::Rotation<double,3> > orientationFunction(orientationFEBasis,
                                                                                                                         x[_1]);
 
   if (dim == dimworld) {
-    CosseratVTKWriter<GridView>::write(gridView,
+    GFE::CosseratVTKWriter<GridView>::write(gridView,
                                        displacementFunction,
                                        orientationFunction,
                                        std::max(LFE_ORDER, GFE_ORDER),
                                        resultPath + "cosserat_homotopy_0_l" + std::to_string(numLevels));
   } else if (dim == 2 && dimworld == 3) {
 #if MIXED_SPACE
-    CosseratVTKWriter<GridView>::write<DeformationFEBasis>(deformationFEBasis, x[_0], resultPath + "cosserat_homotopy_0_l" + std::to_string(numLevels));
+    GFE::CosseratVTKWriter<GridView>::write<DeformationFEBasis>(deformationFEBasis, x[_0], resultPath + "cosserat_homotopy_0_l" + std::to_string(numLevels));
 #else
-    CosseratVTKWriter<GridView>::write<DeformationFEBasis>(deformationFEBasis, x, resultPath + "cosserat_homotopy_0_l" + std::to_string(numLevels));
+    GFE::CosseratVTKWriter<GridView>::write<DeformationFEBasis>(deformationFEBasis, x, resultPath + "cosserat_homotopy_0_l" + std::to_string(numLevels));
 #endif
   }
 
@@ -576,13 +576,13 @@ int main (int argc, char *argv[]) try
     using ScalarDeformationLocalFiniteElement = decltype(compositeBasis.localView().tree().child(_0,0).finiteElement());
     using ScalarRotationLocalFiniteElement = decltype(compositeBasis.localView().tree().child(_1,0).finiteElement());
 
-    using AInterpolationRule = std::tuple<LocalGeodesicFEFunction<dim, double, ScalarDeformationLocalFiniteElement, RealTuple<adouble,3> >,
-        LocalGeodesicFEFunction<dim, double, ScalarRotationLocalFiniteElement, Rotation<adouble,3> > >;
+    using AInterpolationRule = std::tuple<GFE::LocalGeodesicFEFunction<dim, double, ScalarDeformationLocalFiniteElement, GFE::RealTuple<adouble,3> >,
+        GFE::LocalGeodesicFEFunction<dim, double, ScalarRotationLocalFiniteElement, GFE::Rotation<adouble,3> > >;
 
     using ATargetSpace = typename TargetSpace::rebind<adouble>::other;
 
     // The energy on one element
-    auto sumEnergy = std::make_shared<GFE::SumEnergy<CompositeBasis, RealTuple<adouble,3>,Rotation<adouble,3> > >();
+    auto sumEnergy = std::make_shared<GFE::SumEnergy<CompositeBasis, GFE::RealTuple<adouble,3>,GFE::Rotation<adouble,3> > >();
 
     // The actual Cosserat energy
     auto localCosseratEnergy = createCosseratEnergy<CompositeBasis,AInterpolationRule,ATargetSpace,decltype(creator)>(materialParameters,
@@ -591,7 +591,7 @@ int main (int argc, char *argv[]) try
     sumEnergy->addLocalEnergy(localCosseratEnergy);
 
     // The Neumann surface load term
-    auto neumannEnergy = std::make_shared<GFE::NeumannEnergy<CompositeBasis, RealTuple<adouble,3>, Rotation<adouble,3> > >(neumannBoundary,neumannFunction);
+    auto neumannEnergy = std::make_shared<GFE::NeumannEnergy<CompositeBasis, GFE::RealTuple<adouble,3>, GFE::Rotation<adouble,3> > >(neumannBoundary,neumannFunction);
     sumEnergy->addLocalEnergy(neumannEnergy);
 
     // The volume load term
@@ -601,10 +601,10 @@ int main (int argc, char *argv[]) try
     sumEnergy->addLocalEnergy(volumeLoadEnergy);
 
     // The local assembler
-    LocalGeodesicFEADOLCStiffness<CompositeBasis,TargetSpace> localGFEADOLCStiffness(sumEnergy,
+    GFE::LocalGeodesicFEADOLCStiffness<CompositeBasis,TargetSpace> localGFEADOLCStiffness(sumEnergy,
                                                                                      adolcScalarMode);
 
-    MixedGFEAssembler<CompositeBasis,TargetSpace> mixedAssembler(compositeBasis, localGFEADOLCStiffness);
+    GFE::MixedGFEAssembler<CompositeBasis,TargetSpace> mixedAssembler(compositeBasis, localGFEADOLCStiffness);
 
     ////////////////////////////////////////////
     //  Set up the solver
@@ -613,10 +613,10 @@ int main (int argc, char *argv[]) try
 #if MIXED_SPACE
     if (parameterSet.get<std::string>("solvertype", "trustRegion") == "trustRegion")
     {
-      MixedRiemannianTrustRegionSolver<GridType,
+      GFE::MixedRiemannianTrustRegionSolver<GridType,
           CompositeBasis,
-          DeformationFEBasis, RealTuple<double,3>,
-          OrientationFEBasis, Rotation<double,3> > solver;
+          DeformationFEBasis, GFE::RealTuple<double,3>,
+          OrientationFEBasis, GFE::Rotation<double,3> > solver;
       solver.setup(*grid,
                    &mixedAssembler,
                    deformationFEBasis,
@@ -642,7 +642,7 @@ int main (int argc, char *argv[]) try
     else
     {
       //Create BitVector matching the tangential space
-      const int dimRotationTangent = Rotation<double,3>::TangentVector::dimension;
+      const int dimRotationTangent = GFE::Rotation<double,3>::TangentVector::dimension;
       using VectorForBit = MultiTypeBlockVector<std::vector<FieldVector<double,3> >, std::vector<FieldVector<double,dimRotationTangent> > >;
       using BitVector = Solvers::DefaultBitVector_t<VectorForBit>;
       BitVector dirichletDofs;
@@ -658,7 +658,7 @@ int main (int argc, char *argv[]) try
         for (int j = 0; j < dimRotationTangent; j++)
           dirichletDofs[_1][i][j] = orientationDirichletDofs[i][j];
       }
-      GFE::MixedRiemannianProximalNewtonSolver<CompositeBasis, DeformationFEBasis, RealTuple<double,3>, OrientationFEBasis, Rotation<double,3>, BitVector> solver;
+      GFE::MixedRiemannianProximalNewtonSolver<CompositeBasis, DeformationFEBasis, GFE::RealTuple<double,3>, OrientationFEBasis, GFE::Rotation<double,3>, BitVector> solver;
       solver.setup(*grid,
                    &mixedAssembler,
                    x,
@@ -687,10 +687,10 @@ int main (int argc, char *argv[]) try
         dirichletDofsTargetSpace[i][j] = orientationDirichletDofs[i][j-3];
     }
 
-    using GFEAssemblerWrapper = Dune::GFE::GeodesicFEAssemblerWrapper<CompositeBasis, DeformationFEBasis, TargetSpace>;
+    using GFEAssemblerWrapper = GFE::GeodesicFEAssemblerWrapper<CompositeBasis, DeformationFEBasis, TargetSpace>;
     GFEAssemblerWrapper assembler(&mixedAssembler, deformationFEBasis);
     if (parameterSet.get<std::string>("solvertype", "trustRegion") == "trustRegion") {
-      RiemannianTrustRegionSolver<DeformationFEBasis, TargetSpace, GFEAssemblerWrapper> solver;
+      GFE::RiemannianTrustRegionSolver<DeformationFEBasis, TargetSpace, GFEAssemblerWrapper> solver;
       solver.setup(*grid,
                    &assembler,
                    xTargetSpace,
@@ -710,7 +710,7 @@ int main (int argc, char *argv[]) try
       solver.solve();
       xTargetSpace = solver.getSol();
     } else {
-      RiemannianProximalNewtonSolver<DeformationFEBasis, TargetSpace, GFEAssemblerWrapper> solver;
+      GFE::RiemannianProximalNewtonSolver<DeformationFEBasis, TargetSpace, GFEAssemblerWrapper> solver;
       solver.setup(*grid,
                    &assembler,
                    xTargetSpace,
@@ -746,16 +746,16 @@ int main (int argc, char *argv[]) try
     auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,3> >(deformationPowerBasis, displacement);
 
     if (dim == dimworld) {
-      CosseratVTKWriter<GridView>::write(gridView,
+      GFE::CosseratVTKWriter<GridView>::write(gridView,
                                          displacementFunction,
                                          orientationFunction,
                                          std::max(LFE_ORDER, GFE_ORDER),
                                          resultPath + "cosserat_homotopy_" + std::to_string(i+1) + "_l" + std::to_string(numLevels));
     } else if (dim == 2 && dimworld == 3) {
 #if MIXED_SPACE
-      CosseratVTKWriter<GridView>::write<DeformationFEBasis>(deformationFEBasis, x[_0], resultPath + "cosserat_homotopy_" + std::to_string(i+1) + "_l" + std::to_string(numLevels));
+      GFE::CosseratVTKWriter<GridView>::write<DeformationFEBasis>(deformationFEBasis, x[_0], resultPath + "cosserat_homotopy_" + std::to_string(i+1) + "_l" + std::to_string(numLevels));
 #else
-      CosseratVTKWriter<GridView>::write<DeformationFEBasis>(deformationFEBasis, x, resultPath + "cosserat_homotopy_" + std::to_string(i+1) + "_l" + std::to_string(numLevels));
+      GFE::CosseratVTKWriter<GridView>::write<DeformationFEBasis>(deformationFEBasis, x, resultPath + "cosserat_homotopy_" + std::to_string(i+1) + "_l" + std::to_string(numLevels));
 #endif
     }
   }
diff --git a/src/cosserat-rod.cc b/src/cosserat-rod.cc
index 27daac7f..86d8a1f1 100644
--- a/src/cosserat-rod.cc
+++ b/src/cosserat-rod.cc
@@ -42,7 +42,7 @@
 using namespace Dune;
 using namespace Dune::Indices;
 
-using TargetSpace = GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> >;
+using TargetSpace = GFE::ProductManifold<GFE::RealTuple<double,3>,GFE::Rotation<double,3> >;
 
 const int blocksize = TargetSpace::TangentVector::dimension;
 
@@ -220,7 +220,7 @@ int main (int argc, char *argv[]) try
   //////////////////////////////////////////////
 
   using ATargetSpace = TargetSpace::rebind<adouble>::other;
-  using GeodesicInterpolationRule  = LocalGeodesicFEFunction<1, double, ScalarBasis::LocalView::Tree::FiniteElement, ATargetSpace>;
+  using GeodesicInterpolationRule  = GFE::LocalGeodesicFEFunction<1, double, ScalarBasis::LocalView::Tree::FiniteElement, ATargetSpace>;
   using ProjectedInterpolationRule = GFE::LocalProjectedFEFunction<1, double, ScalarBasis::LocalView::Tree::FiniteElement, ATargetSpace>;
 
   // Assembler using ADOL-C
@@ -243,16 +243,16 @@ int main (int argc, char *argv[]) try
   else
     DUNE_THROW(Exception, "Unknown interpolation method " << parameterSet["interpolationMethod"] << " requested!");
 
-  LocalGeodesicFEADOLCStiffness<ScalarBasis,
+  GFE::LocalGeodesicFEADOLCStiffness<ScalarBasis,
       TargetSpace> localStiffness(localRodEnergy);
 
-  GeodesicFEAssembler<ScalarBasis,TargetSpace> rodAssembler(gridView, localStiffness);
+  GFE::GeodesicFEAssembler<ScalarBasis,TargetSpace> rodAssembler(gridView, localStiffness);
 
   /////////////////////////////////////////////
   //   Create a solver for the rod problem
   /////////////////////////////////////////////
 
-  RiemannianTrustRegionSolver<ScalarBasis,TargetSpace> rodSolver;
+  GFE::RiemannianTrustRegionSolver<ScalarBasis,TargetSpace> rodSolver;
 
   rodSolver.setup(grid,
                   &rodAssembler,
@@ -309,16 +309,16 @@ int main (int argc, char *argv[]) try
   auto displacementFunction = Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,3> >(worldBasis, displacement);
 
   // Copy the orientation part of the configuration; the CosseratVTKWriter wants it that way
-  std::vector<Rotation<double,3> > orientationConfiguration(x.size());
+  std::vector<GFE::Rotation<double,3> > orientationConfiguration(x.size());
   for (size_t i=0; i<x.size(); ++i)
     orientationConfiguration[i] = x[i][_1];
 
-  using RotationInterpolationRule  = LocalGeodesicFEFunction<1, double, ScalarBasis::LocalView::Tree::FiniteElement, Rotation<double,3> >;
+  using RotationInterpolationRule = GFE::LocalGeodesicFEFunction<1, double, ScalarBasis::LocalView::Tree::FiniteElement, GFE::Rotation<double,3> >;
 
-  GFE::EmbeddedGlobalGFEFunction<ScalarBasis, RotationInterpolationRule,Rotation<double,3> > orientationFunction(scalarBasis,
+  GFE::EmbeddedGlobalGFEFunction<ScalarBasis, RotationInterpolationRule,GFE::Rotation<double,3> > orientationFunction(scalarBasis,
                                                                                                                  orientationConfiguration);
 
-  CosseratVTKWriter<GridView>::write(gridView,
+  GFE::CosseratVTKWriter<GridView>::write(gridView,
                                      displacementFunction,
                                      orientationFunction,
                                      order,
diff --git a/src/film-on-substrate-stress-plot.cc b/src/film-on-substrate-stress-plot.cc
index 29524617..7ccbee58 100644
--- a/src/film-on-substrate-stress-plot.cc
+++ b/src/film-on-substrate-stress-plot.cc
@@ -176,7 +176,7 @@ int main (int argc, char *argv[]) try
   auto deformationMap = Dune::GFE::transformFileToMap<dim>(pathToOutput + parameterSet.get<std::string>("deformationOutput"));
   std::cout << "... done: The basis has " << basisOrderD.size() << " elements and the defomation file has " << deformationMap.size() << " entries." << std::endl;
 
-  const auto dimRotation = Rotation<double,dim>::embeddedDim;
+  const auto dimRotation = GFE::Rotation<double,dim>::embeddedDim;
   std::unordered_map<std::string, FieldVector<double,dimRotation> > rotationMap;
   if (parameterSet.hasKey("rotationOutput")) {
     std::cout << "Reading in rotation file ("  << "order is "  << rotationOrder  << "): " << pathToOutput + parameterSet.get<std::string>("rotationOutput") << std::endl;
@@ -222,7 +222,7 @@ int main (int argc, char *argv[]) try
     }
   }
 
-  using RotationVector = std::vector<Rotation<double,dim> >;
+  using RotationVector = std::vector<GFE::Rotation<double,dim> >;
   RotationVector rot;
   rot.resize(basisOrderR.size());
   DisplacementVector xOrderR;
@@ -240,7 +240,7 @@ int main (int argc, char *argv[]) try
   for (std::size_t i = 0; i < basisOrderR.size(); i++) {
     std::stringstream stream;
     stream << xOrderR[i];
-    Rotation<double,dim> rotation(rotationMap.at(stream.str()));
+    GFE::Rotation<double,dim> rotation(rotationMap.at(stream.str()));
     FieldMatrix<double,dim,dim> rotationMatrix(0);
     rotation.matrix(rotationMatrix);
     rot[i].set(rotationMatrix);
@@ -253,7 +253,7 @@ int main (int argc, char *argv[]) try
   /////////////////////////////////////////////////////////////
   int quadOrder = parameterSet.hasKey("quadOrder") ? parameterSet.get<int>("quadOrder") : 4;
 
-  auto stressAssembler = GFE::SurfaceCosseratStressAssembler<decltype(basisOrderD),decltype(basisOrderR), FieldVector<double,dim>, Rotation<double,dim> >
+  auto stressAssembler = GFE::SurfaceCosseratStressAssembler<decltype(basisOrderD),decltype(basisOrderR), FieldVector<double,dim>, GFE::Rotation<double,dim> >
                            (basisOrderD, basisOrderR);
 
 
diff --git a/src/film-on-substrate.cc b/src/film-on-substrate.cc
index a9633e60..8cf8189d 100644
--- a/src/film-on-substrate.cc
+++ b/src/film-on-substrate.cc
@@ -231,9 +231,9 @@ int main (int argc, char *argv[]) try
 
   using namespace Dune::Indices;
 
-  typedef std::vector<RealTuple<double,dim> > DisplacementVector;
-  typedef std::vector<Rotation<double,dim> > RotationVector;
-  const int dimRotation = Rotation<double,dim>::embeddedDim;
+  typedef std::vector<GFE::RealTuple<double,dim> > DisplacementVector;
+  typedef std::vector<GFE::Rotation<double,dim> > RotationVector;
+  const int dimRotation = GFE::Rotation<double,dim>::embeddedDim;
   typedef TupleVector<DisplacementVector, RotationVector> SolutionType;
 
   /////////////////////////////////////////////////////////////
@@ -328,7 +328,7 @@ int main (int argc, char *argv[]) try
 #endif
 
   //Create BitVector matching the tangential space
-  const int dimRotationTangent = Rotation<double,dim>::TangentVector::dimension;
+  const int dimRotationTangent = GFE::Rotation<double,dim>::TangentVector::dimension;
   typedef MultiTypeBlockVector<std::vector<FieldVector<double,dim> >, std::vector<FieldVector<double,dimRotationTangent> > > VectorForBit;
   typedef Solvers::DefaultBitVector_t<VectorForBit> BitVector;
 
@@ -515,7 +515,7 @@ int main (int argc, char *argv[]) try
     if(!elasticDensity)
       DUNE_THROW(Exception, "Error: Selected energy not available!");
 
-    using ActiveRigidBodyMotion = GFE::ProductManifold<RealTuple<adouble,dim>, Rotation<adouble,dim> >;
+    using ActiveRigidBodyMotion = GFE::ProductManifold<GFE::RealTuple<adouble,dim>, GFE::Rotation<adouble,dim> >;
 
     // Wrap dune-elasticity density as dune-gfe density
     using Element = GridView::Codim<0>::Entity;
@@ -523,13 +523,13 @@ int main (int argc, char *argv[]) try
 
 
     // Select which type of geometric interpolation to use
-    using LocalDeformationInterpolationRule = LocalGeodesicFEFunction<dim, GridType::ctype, decltype(deformationFEBasis.localView().tree().finiteElement()), RealTuple<adouble,dim> >;
-    using LocalOrientationInterpolationRule = LocalGeodesicFEFunction<dim, GridType::ctype, decltype(orientationFEBasis.localView().tree().finiteElement()), Rotation<adouble,dim> >;
+    using LocalDeformationInterpolationRule = GFE::LocalGeodesicFEFunction<dim, GridType::ctype, decltype(deformationFEBasis.localView().tree().finiteElement()), GFE::RealTuple<adouble,dim> >;
+    using LocalOrientationInterpolationRule = GFE::LocalGeodesicFEFunction<dim, GridType::ctype, decltype(orientationFEBasis.localView().tree().finiteElement()), GFE::Rotation<adouble,dim> >;
 
     using LocalInterpolationRule = std::tuple<LocalDeformationInterpolationRule,LocalOrientationInterpolationRule>;
 
     auto elasticEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, LocalInterpolationRule, ActiveRigidBodyMotion> >(elasticDensityWrapped);
-    auto neumannEnergy = std::make_shared<GFE::NeumannEnergy<CompositeBasis, RealTuple<ValueType,targetDim>, Rotation<ValueType,dim> > >(neumannBoundary,neumannFunctionPtr);
+    auto neumannEnergy = std::make_shared<GFE::NeumannEnergy<CompositeBasis, GFE::RealTuple<ValueType,targetDim>, GFE::Rotation<ValueType,dim> > >(neumannBoundary,neumannFunctionPtr);
 
     // The energy of the surface shell
     using Intersection = typename GridView::Intersection;
@@ -544,27 +544,27 @@ int main (int argc, char *argv[]) try
       &surfaceShellBoundary,
       stressFreeShellFunction);
 
-    using RBM = GFE::ProductManifold<RealTuple<double, dim>,Rotation<double,dim> >;
+    using RBM = GFE::ProductManifold<GFE::RealTuple<double, dim>,GFE::Rotation<double,dim> >;
 
-    auto sumEnergy = std::make_shared<GFE::SumEnergy<CompositeBasis, RealTuple<ValueType,targetDim>, Rotation<ValueType,targetDim> > >();
+    auto sumEnergy = std::make_shared<GFE::SumEnergy<CompositeBasis, GFE::RealTuple<ValueType,targetDim>, GFE::Rotation<ValueType,targetDim> > >();
     sumEnergy->addLocalEnergy(neumannEnergy);
     sumEnergy->addLocalEnergy(elasticEnergy);
     sumEnergy->addLocalEnergy(surfaceCosseratEnergy);
 
-    LocalGeodesicFEADOLCStiffness<CompositeBasis,RBM> localGFEADOLCStiffness(sumEnergy);
-    MixedGFEAssembler<CompositeBasis,RBM> mixedAssembler(compositeBasis, localGFEADOLCStiffness);
+    GFE::LocalGeodesicFEADOLCStiffness<CompositeBasis,RBM> localGFEADOLCStiffness(sumEnergy);
+    GFE::MixedGFEAssembler<CompositeBasis,RBM> mixedAssembler(compositeBasis, localGFEADOLCStiffness);
 
     ////////////////////////////////////////////////////////
     //   Set Dirichlet values
     ////////////////////////////////////////////////////////
 
-    BitSetVector<RealTuple<double,dim>::TangentVector::dimension> deformationDirichletDofs(dirichletDofs[_0].size(), false);
+    BitSetVector<GFE::RealTuple<double,dim>::TangentVector::dimension> deformationDirichletDofs(dirichletDofs[_0].size(), false);
     for (std::size_t i = 0; i < dirichletDofs[_0].size(); i++)
-      for (int j = 0; j < RealTuple<double,dim>::TangentVector::dimension; j++)
+      for (int j = 0; j < GFE::RealTuple<double,dim>::TangentVector::dimension; j++)
         deformationDirichletDofs[i][j] = dirichletDofs[_0][i][j];
-    BitSetVector<Rotation<double,dim>::TangentVector::dimension> orientationDirichletDofs(dirichletDofs[_1].size(), false);
+    BitSetVector<GFE::Rotation<double,dim>::TangentVector::dimension> orientationDirichletDofs(dirichletDofs[_1].size(), false);
     for (std::size_t i = 0; i < dirichletDofs[_1].size(); i++)
-      for (int j = 0; j < Rotation<double,dim>::TangentVector::dimension; j++)
+      for (int j = 0; j < GFE::Rotation<double,dim>::TangentVector::dimension; j++)
         orientationDirichletDofs[i][j] = dirichletDofs[_1][i][j];
 
     Python::Reference dirichletValuesClass = Python::import(parameterSet.get<std::string>("dirichletValues"));
@@ -619,7 +619,7 @@ int main (int argc, char *argv[]) try
 
     if (parameterSet.get<std::string>("solvertype", "trustRegion") == "trustRegion") {
 #if MIXED_SPACE
-      MixedRiemannianTrustRegionSolver<GridType, CompositeBasis, DeformationFEBasis, RealTuple<double,dim>, OrientationFEBasis, Rotation<double,dim> > solver;
+      GFE::MixedRiemannianTrustRegionSolver<GridType, CompositeBasis, DeformationFEBasis, GFE::RealTuple<double,dim>, OrientationFEBasis, GFE::Rotation<double,dim> > solver;
       solver.setup(*grid,
                    &mixedAssembler,
                    deformationFEBasis,
@@ -642,7 +642,7 @@ int main (int argc, char *argv[]) try
       solver.solve();
       x = solver.getSol();
 #else
-      RiemannianTrustRegionSolver<DeformationFEBasis, RBM, GFEAssemblerWrapper> solver;
+      GFE::RiemannianTrustRegionSolver<DeformationFEBasis, RBM, GFEAssemblerWrapper> solver;
       solver.setup(*grid,
                    &assembler,
                    xRBM,
@@ -669,7 +669,7 @@ int main (int argc, char *argv[]) try
     } else { //parameterSet.get<std::string>("solvertype") == "proximalNewton"
 
 #if MIXED_SPACE
-      GFE::MixedRiemannianProximalNewtonSolver<CompositeBasis, DeformationFEBasis, RealTuple<double,dim>, OrientationFEBasis, Rotation<double,dim>, BitVector> solver;
+      GFE::MixedRiemannianProximalNewtonSolver<CompositeBasis, DeformationFEBasis, GFE::RealTuple<double,dim>, OrientationFEBasis, GFE::Rotation<double,dim>, BitVector> solver;
       solver.setup(*grid,
                    &mixedAssembler,
                    x,
@@ -682,7 +682,7 @@ int main (int argc, char *argv[]) try
       solver.solve();
       x = solver.getSol();
 #else
-      RiemannianProximalNewtonSolver<DeformationFEBasis, RBM, GFEAssemblerWrapper> solver;
+      GFE::RiemannianProximalNewtonSolver<DeformationFEBasis, RBM, GFEAssemblerWrapper> solver;
       solver.setup(*grid,
                    &assembler,
                    xRBM,
diff --git a/src/gradient-flow.cc b/src/gradient-flow.cc
index a4197e7d..bb88f764 100644
--- a/src/gradient-flow.cc
+++ b/src/gradient-flow.cc
@@ -45,17 +45,19 @@
 #include <dune/gfe/densities/harmonicdensity.hh>
 #include <dune/gfe/spaces/unitvector.hh>
 
+using namespace Dune;
+
 // grid dimension
 const int dim = 1;
 const int dimworld = dim;
 
 // Image space of the geodesic fe functions
-// typedef Rotation<double,2> TargetSpace;
-// typedef Rotation<double,3> TargetSpace;
-// typedef UnitVector<double,2> TargetSpace;
-typedef UnitVector<double,3> TargetSpace;
-// typedef UnitVector<double,4> TargetSpace;
-// typedef RealTuple<double,1> TargetSpace;
+// typedef GFE::Rotation<double,2> TargetSpace;
+// typedef GFE::Rotation<double,3> TargetSpace;
+// typedef GFE::UnitVector<double,2> TargetSpace;
+typedef GFE::UnitVector<double,3> TargetSpace;
+// typedef GFE::UnitVector<double,4> TargetSpace;
+// typedef GFE::RealTuple<double,1> TargetSpace;
 
 // Tangent vector of the image space
 const int blocksize = TargetSpace::TangentVector::dimension;
@@ -63,7 +65,6 @@ const int blocksize = TargetSpace::TangentVector::dimension;
 // Approximation order of the finite element space
 const int order = 1;
 
-using namespace Dune;
 
 
 int main (int argc, char *argv[]) try
@@ -205,10 +206,10 @@ int main (int argc, char *argv[]) try
   // Assembler using ADOL-C
   typedef TargetSpace::rebind<adouble>::other ATargetSpace;
 
-  auto l2DistanceSquaredEnergy = std::make_shared<L2DistanceSquaredEnergy<FEBasis, ATargetSpace> >();
+  auto l2DistanceSquaredEnergy = std::make_shared<GFE::L2DistanceSquaredEnergy<FEBasis, ATargetSpace> >();
 
   std::vector<std::shared_ptr<GFE::LocalEnergy<FEBasis,ATargetSpace> > > addends(2);
-  using GeodesicInterpolationRule  = LocalGeodesicFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, ATargetSpace>;
+  using GeodesicInterpolationRule = GFE::LocalGeodesicFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, ATargetSpace>;
 
   auto harmonicDensity = std::make_shared<GFE::HarmonicDensity<GridType::Codim<0>::Entity, ATargetSpace> >();
   addends[0] = std::make_shared<GFE::LocalIntegralEnergy<FEBasis, GeodesicInterpolationRule, ATargetSpace> >(harmonicDensity);
@@ -216,17 +217,17 @@ int main (int argc, char *argv[]) try
 
   std::vector<double> weights = {1.0, 1.0/(2*timeStepSize)};
 
-  auto sumEnergy = std::make_shared< WeightedSumEnergy<FEBasis, ATargetSpace> >(addends, weights);
+  auto sumEnergy = std::make_shared< GFE::WeightedSumEnergy<FEBasis, ATargetSpace> >(addends, weights);
 
-  LocalGeodesicFEADOLCStiffness<FEBasis,TargetSpace> localGFEADOLCStiffness(sumEnergy);
+  GFE::LocalGeodesicFEADOLCStiffness<FEBasis,TargetSpace> localGFEADOLCStiffness(sumEnergy);
 
-  GeodesicFEAssembler<FEBasis,TargetSpace> assembler(feBasis, localGFEADOLCStiffness);
+  GFE::GeodesicFEAssembler<FEBasis,TargetSpace> assembler(feBasis, localGFEADOLCStiffness);
 
   ///////////////////////////////////////////////////
   //   Create a Riemannian trust-region solver
   ///////////////////////////////////////////////////
 
-  RiemannianTrustRegionSolver<FEBasis,TargetSpace> solver;
+  GFE::RiemannianTrustRegionSolver<FEBasis,TargetSpace> solver;
   solver.setup(*grid,
                &assembler,
                x,
diff --git a/src/harmonicmaps.cc b/src/harmonicmaps.cc
index 2eb356ba..5f198e9e 100644
--- a/src/harmonicmaps.cc
+++ b/src/harmonicmaps.cc
@@ -57,20 +57,21 @@
 const int dim = 2;
 const int dimworld = 2;
 
+using namespace Dune;
+
 // Image space of the geodesic fe functions
-// typedef Rotation<double,2> TargetSpace;
-// typedef Rotation<double,3> TargetSpace;
-// typedef UnitVector<double,2> TargetSpace;
-typedef UnitVector<double,3> TargetSpace;
-// typedef UnitVector<double,4> TargetSpace;
-// typedef RealTuple<double,1> TargetSpace;
+// typedef GFE::Rotation<double,2> TargetSpace;
+// typedef GFE::Rotation<double,3> TargetSpace;
+// typedef GFE::UnitVector<double,2> TargetSpace;
+typedef GFE::UnitVector<double,3> TargetSpace;
+// typedef GFE::UnitVector<double,4> TargetSpace;
+// typedef GFE::RealTuple<double,1> TargetSpace;
 
 // Tangent vector of the image space
 const int blocksize = TargetSpace::TangentVector::dimension;
 
 const int order = 1;
 
-using namespace Dune;
 
 template <typename Writer, typename Basis, typename SolutionType>
 void fillVTKWriter(Writer& vtkWriter, const Basis& feBasis, const SolutionType& x, std::string filename)
@@ -80,7 +81,7 @@ void fillVTKWriter(Writer& vtkWriter, const Basis& feBasis, const SolutionType&
   for (size_t i=0; i<x.size(); i++)
     xEmbedded[i] = x[i].globalCoordinates();
 
-  if constexpr (std::is_same<TargetSpace, Rotation<double,3> >::value)
+  if constexpr (std::is_same<TargetSpace, GFE::Rotation<double,3> >::value)
   {
     std::array<BlockVector<FieldVector<double,3> >,3> director;
     for (int i=0; i<3; i++)
@@ -287,7 +288,7 @@ int main (int argc, char *argv[])
     DUNE_THROW(Exception, "Unknown energy type '" << energy << "'");
 
   // Next: The local energy, i.e., the integral of the density over one element
-  using GeodesicInterpolationRule  = LocalGeodesicFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, ATargetSpace>;
+  using GeodesicInterpolationRule  = GFE::LocalGeodesicFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, ATargetSpace>;
   using ProjectedInterpolationRule = GFE::LocalProjectedFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, ATargetSpace>;
 
   std::shared_ptr<GFE::LocalEnergy<FEBasis,ATargetSpace> > localEnergy;
@@ -300,15 +301,15 @@ int main (int argc, char *argv[])
     DUNE_THROW(Exception, "Unknown interpolation method " << parameterSet["interpolationMethod"] << " requested!");
 
   // Compute local tangent problems by applying ADOL-C directly to the energy on the element
-  LocalGeodesicFEADOLCStiffness<FEBasis,TargetSpace> localGFEADOLCStiffness(localEnergy);
+  GFE::LocalGeodesicFEADOLCStiffness<FEBasis,TargetSpace> localGFEADOLCStiffness(localEnergy);
 
-  GeodesicFEAssembler<FEBasis,TargetSpace> assembler(feBasis, localGFEADOLCStiffness);
+  GFE::GeodesicFEAssembler<FEBasis,TargetSpace> assembler(feBasis, localGFEADOLCStiffness);
 
   // /////////////////////////////////////////////////
   //   Create a Riemannian trust-region solver
   // /////////////////////////////////////////////////
 
-  RiemannianTrustRegionSolver<FEBasis,TargetSpace> solver;
+  GFE::RiemannianTrustRegionSolver<FEBasis,TargetSpace> solver;
   solver.setup(*grid,
                &assembler,
                x,
diff --git a/src/simofoxshell.cc b/src/simofoxshell.cc
index 968ef0b4..3c47af5c 100644
--- a/src/simofoxshell.cc
+++ b/src/simofoxshell.cc
@@ -75,7 +75,7 @@ int main(int argc, char *argv[]) try
     << std::endl;
 
   using namespace Dune::Indices;
-  using SolutionType = TupleVector<std::vector<RealTuple<double,3> >, std::vector<UnitVector<double,3> > >;
+  using SolutionType = TupleVector<std::vector<GFE::RealTuple<double,3> >, std::vector<GFE::UnitVector<double,3> > >;
 
   // parse data file
   ParameterTree parameterSet;
@@ -325,12 +325,12 @@ int main(int argc, char *argv[]) try
                     neumannFunction,
                     nullptr, x0);
 
-    using TargetSpace = Dune::GFE::ProductManifold<RealTuple<double,3>,UnitVector<double,3> >;
+    using TargetSpace = GFE::ProductManifold<GFE::RealTuple<double,3>,GFE::UnitVector<double,3> >;
 
-    LocalGeodesicFEADOLCStiffness<decltype(compositeBasis),
+    GFE::LocalGeodesicFEADOLCStiffness<decltype(compositeBasis),
         TargetSpace> localGFEADOLCStiffness(simoFoxEnergyLocalStiffness);
 
-    MixedGFEAssembler<decltype(compositeBasis),TargetSpace> assembler(compositeBasis, localGFEADOLCStiffness);
+    GFE::MixedGFEAssembler<decltype(compositeBasis),TargetSpace> assembler(compositeBasis, localGFEADOLCStiffness);
     ////////////////////////////////////////////////////////
     //   Set Dirichlet values
     ////////////////////////////////////////////////////////
@@ -359,10 +359,10 @@ int main(int argc, char *argv[]) try
     // /////////////////////////////////////////////////
     if (parameterSet.get<std::string>("solvertype", "trustRegion") == "trustRegion") {
 
-      MixedRiemannianTrustRegionSolver<Grid,
+      GFE::MixedRiemannianTrustRegionSolver<Grid,
           decltype(compositeBasis),
-          MidsurfaceFEBasis, RealTuple<double,3>,
-          DirectorFEBasis, UnitVector<double,3> > solver;
+          MidsurfaceFEBasis, GFE::RealTuple<double,3>,
+          DirectorFEBasis, GFE::UnitVector<double,3> > solver;
 
       solver.setup(*grid,
                    &assembler,
@@ -403,9 +403,9 @@ int main(int argc, char *argv[]) try
         for (int j = 3; j < TargetSpace::TangentVector::dimension; j ++)
           dirichletDofsTargetSpace[i][j] = orientationDirichletDofs[i][j-3];
       }
-      using GFEAssemblerWrapper = Dune::GFE::GeodesicFEAssemblerWrapper<decltype(compositeBasis), MidsurfaceFEBasis, TargetSpace>;
+      using GFEAssemblerWrapper = GFE::GeodesicFEAssemblerWrapper<decltype(compositeBasis), MidsurfaceFEBasis, TargetSpace>;
       GFEAssemblerWrapper assemblerNotMixed(&assembler, midsurfaceFEBasis);
-      RiemannianProximalNewtonSolver<MidsurfaceFEBasis, TargetSpace, GFEAssemblerWrapper> solver;
+      GFE::RiemannianProximalNewtonSolver<MidsurfaceFEBasis, TargetSpace, GFEAssemblerWrapper> solver;
       solver.setup(*grid,
                    &assemblerNotMixed,
                    xTargetSpace,
diff --git a/test/averagedistanceassemblertest.cc b/test/averagedistanceassemblertest.cc
index 144e9b29..e0bebced 100644
--- a/test/averagedistanceassemblertest.cc
+++ b/test/averagedistanceassemblertest.cc
@@ -104,7 +104,7 @@ void testPoint(const std::vector<TargetSpace>& corners,
                const TargetSpace& argument)
 {
   // create the assembler
-  AverageDistanceAssembler<TargetSpace> assembler(corners, weights);
+  GFE::AverageDistanceAssembler<TargetSpace> assembler(corners, weights);
 
   // test the functional
   double value = assembler.value(argument);
@@ -180,7 +180,7 @@ void testWeightSet(const std::vector<TargetSpace>& corners,
 
 void testRealTuples()
 {
-  typedef RealTuple<double,1> TargetSpace;
+  typedef GFE::RealTuple<double,1> TargetSpace;
 
   std::vector<TargetSpace> corners = {TargetSpace(1),
                                       TargetSpace(2),
@@ -196,7 +196,7 @@ void testRealTuples()
 
 void testUnitVectors()
 {
-  typedef UnitVector<double,3> TargetSpace;
+  typedef GFE::UnitVector<double,3> TargetSpace;
 
   std::vector<TargetSpace> corners(dim+1);
 
@@ -214,12 +214,12 @@ void testUnitVectors()
 
 void testRotations()
 {
-  typedef Rotation<double,3> TargetSpace;
+  typedef GFE::Rotation<double,3> TargetSpace;
 
   std::vector<TargetSpace> corners(dim+1);
-  corners[0] = Rotation<double,3>({1,0,0}, 0.1);
-  corners[1] = Rotation<double,3>({0,1,0}, 0.1);
-  corners[2] = Rotation<double,3>({0,0,1}, 0.1);
+  corners[0] = GFE::Rotation<double,3>({1,0,0}, 0.1);
+  corners[1] = GFE::Rotation<double,3>({0,1,0}, 0.1);
+  corners[2] = GFE::Rotation<double,3>({0,0,1}, 0.1);
 
   TargetSpace argument = corners[0];
   testWeightSet(corners, argument);
@@ -232,12 +232,12 @@ void testRotations()
 
 void testProductManifold()
 {
-  typedef Dune::GFE::ProductManifold<RealTuple<double,5>,UnitVector<double,3>, Rotation<double,3> > TargetSpace;
+  typedef GFE::ProductManifold<GFE::RealTuple<double,5>,GFE::UnitVector<double,3>, GFE::Rotation<double,3> > TargetSpace;
 
   std::vector<TargetSpace> corners(dim+1);
 
   std::generate(corners.begin(), corners.end(), []()  {
-    return Dune::GFE::randomFieldVector<typename TargetSpace::field_type,TargetSpace::CoordinateType::dimension>(0.9,1.1);
+    return GFE::randomFieldVector<typename TargetSpace::field_type,TargetSpace::CoordinateType::dimension>(0.9,1.1);
   });
 
   TargetSpace argument = corners[0];
diff --git a/test/cosseratcontinuumtest.cc b/test/cosseratcontinuumtest.cc
index f0049e73..a17db46e 100644
--- a/test/cosseratcontinuumtest.cc
+++ b/test/cosseratcontinuumtest.cc
@@ -45,7 +45,7 @@ int main (int argc, char *argv[])
 
   MPIHelper::instance(argc, argv);
 
-  using SolutionType = TupleVector<std::vector<RealTuple<double,dim> >, std::vector<Rotation<double,dim> > >;
+  using SolutionType = TupleVector<std::vector<GFE::RealTuple<double,dim> >, std::vector<GFE::Rotation<double,dim> > >;
 
   // solver settings
   const double tolerance                = 1e-6;
@@ -84,7 +84,7 @@ int main (int argc, char *argv[])
   using namespace Dune::Indices;
   using namespace Functions::BasisFactory;
 
-  const int dimRotation = Rotation<double,dim>::embeddedDim;
+  const int dimRotation = GFE::Rotation<double,dim>::embeddedDim;
   auto compositeBasis = makeBasis(
     gridView,
     composite(
@@ -198,15 +198,15 @@ int main (int argc, char *argv[])
   //  Create an assembler
   // ////////////////////////////
 
-  using RigidBodyMotion = GFE::ProductManifold<RealTuple<double,dim>,Rotation<double,dim> >;
-  using ARigidBodyMotion = GFE::ProductManifold<RealTuple<adouble,dim>,Rotation<adouble,dim> >;
+  using RigidBodyMotion = GFE::ProductManifold<GFE::RealTuple<double,dim>,GFE::Rotation<double,dim> >;
+  using ARigidBodyMotion = GFE::ProductManifold<GFE::RealTuple<adouble,dim>,GFE::Rotation<adouble,dim> >;
 
-  auto sumEnergy = std::make_shared<GFE::SumEnergy<CompositeBasis, RealTuple<adouble,dim>,Rotation<adouble,dim> > >();
-  auto neumannEnergy = std::make_shared<GFE::NeumannEnergy<CompositeBasis, RealTuple<adouble,dim>, Rotation<adouble,dim> > >(neumannBoundary,neumannFunction);
+  auto sumEnergy = std::make_shared<GFE::SumEnergy<CompositeBasis, GFE::RealTuple<adouble,dim>,GFE::Rotation<adouble,dim> > >();
+  auto neumannEnergy = std::make_shared<GFE::NeumannEnergy<CompositeBasis, GFE::RealTuple<adouble,dim>, GFE::Rotation<adouble,dim> > >(neumannBoundary,neumannFunction);
 
   // Select which type of geometric interpolation to use
-  using LocalDeformationInterpolationRule = LocalGeodesicFEFunction<dim, GridType::ctype, decltype(deformationFEBasis.localView().tree().finiteElement()), RealTuple<adouble,dim> >;
-  using LocalOrientationInterpolationRule = LocalGeodesicFEFunction<dim, GridType::ctype, decltype(orientationFEBasis.localView().tree().finiteElement()), Rotation<adouble,dim> >;
+  using LocalDeformationInterpolationRule = GFE::LocalGeodesicFEFunction<dim, GridType::ctype, decltype(deformationFEBasis.localView().tree().finiteElement()), GFE::RealTuple<adouble,dim> >;
+  using LocalOrientationInterpolationRule = GFE::LocalGeodesicFEFunction<dim, GridType::ctype, decltype(orientationFEBasis.localView().tree().finiteElement()), GFE::Rotation<adouble,dim> >;
 
   using LocalInterpolationRule = std::tuple<LocalDeformationInterpolationRule,LocalOrientationInterpolationRule>;
 
@@ -216,13 +216,13 @@ int main (int argc, char *argv[])
   sumEnergy->addLocalEnergy(bulkCosseratEnergy);
   sumEnergy->addLocalEnergy(neumannEnergy);
 
-  LocalGeodesicFEADOLCStiffness<CompositeBasis,RigidBodyMotion> localGFEADOLCStiffness(sumEnergy);
-  MixedGFEAssembler<CompositeBasis,RigidBodyMotion> mixedAssembler(compositeBasis, localGFEADOLCStiffness);
+  GFE::LocalGeodesicFEADOLCStiffness<CompositeBasis,RigidBodyMotion> localGFEADOLCStiffness(sumEnergy);
+  GFE::MixedGFEAssembler<CompositeBasis,RigidBodyMotion> mixedAssembler(compositeBasis, localGFEADOLCStiffness);
 
-  MixedRiemannianTrustRegionSolver<GridType,
+  GFE::MixedRiemannianTrustRegionSolver<GridType,
       CompositeBasis,
-      DeformationFEBasis, RealTuple<double,dim>,
-      OrientationFEBasis, Rotation<double,dim> > solver;
+      DeformationFEBasis, GFE::RealTuple<double,dim>,
+      OrientationFEBasis, GFE::Rotation<double,dim> > solver;
   solver.setup(*grid,
                &mixedAssembler,
                deformationFEBasis,
diff --git a/test/cosseratrodenergytest.cc b/test/cosseratrodenergytest.cc
index 0bdd3b69..8a9c6e0f 100644
--- a/test/cosseratrodenergytest.cc
+++ b/test/cosseratrodenergytest.cc
@@ -18,7 +18,7 @@ using namespace Dune::Indices;
 int main (int argc, char *argv[]) try
 {
   // Type used for algebraic rod configurations
-  using TargetSpace = GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> >;
+  using TargetSpace = GFE::ProductManifold<GFE::RealTuple<double,3>,GFE::Rotation<double,3> >;
   using SolutionType = std::vector<TargetSpace>;
 
   // Problem settings
@@ -45,12 +45,12 @@ int main (int argc, char *argv[]) try
   {
     double s = double(i)/(x.size()-1);
     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)) );
+    x[i][_1] = GFE::Rotation<double,3>::identity();
+    //x[i].q = GFE::Quaternion<double>(zAxis, (double(i)*M_PI)/(2*(x.size()-1)) );
   }
 
   FieldVector<double,3> zAxis(0);  zAxis[2]=1;
-  x.back()[_1] = Rotation<double,3>(zAxis, M_PI/4);
+  x.back()[_1] = GFE::Rotation<double,3>(zAxis, M_PI/4);
 
   // /////////////////////////////////////////////////////////////////////
   //   Create a second, rotated copy of the configuration
@@ -59,7 +59,7 @@ int main (int argc, char *argv[]) try
   FieldVector<double,3> displacement {0, 1, 0};
 
   FieldVector<double,3> axis = {1,0,0};
-  Rotation<double,3> rotation(axis,M_PI/2);
+  GFE::Rotation<double,3> rotation(axis,M_PI/2);
 
   SolutionType rotatedX = x;
 
@@ -71,7 +71,7 @@ int main (int argc, char *argv[]) try
     rotatedX[i][_1] = rotation.mult(x[i][_1]);
   }
 
-  using GeodesicInterpolationRule  = LocalGeodesicFEFunction<1, double,
+  using GeodesicInterpolationRule  = GFE::LocalGeodesicFEFunction<1, double,
       FEBasis::LocalView::Tree::FiniteElement,
       TargetSpace>;
 
@@ -87,7 +87,7 @@ int main (int argc, char *argv[]) try
     auto idx = gridView.indexSet().index(vertex);
 
     referenceConfiguration[idx][_0] = {0.0, 0.0, vertex.geometry().corner(0)[0]};
-    referenceConfiguration[idx][_1] = Rotation<double,3>::identity();
+    referenceConfiguration[idx][_1] = GFE::Rotation<double,3>::identity();
   }
 
   localRodEnergy.setReferenceConfiguration(referenceConfiguration);
diff --git a/test/cosseratrodtest.cc b/test/cosseratrodtest.cc
index 037b02aa..c3758b8a 100644
--- a/test/cosseratrodtest.cc
+++ b/test/cosseratrodtest.cc
@@ -30,7 +30,7 @@
 using namespace Dune;
 using namespace Dune::Indices;
 
-using TargetSpace = GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> >;
+using TargetSpace = GFE::ProductManifold<GFE::RealTuple<double,3>,GFE::Rotation<double,3> >;
 
 const int blocksize = TargetSpace::TangentVector::dimension;
 
@@ -80,7 +80,7 @@ int main (int argc, char *argv[]) try
   for (std::size_t i=0; i<referenceConfiguration.size(); i++)
   {
     referenceConfiguration[i][_0] = {0.0, 0.0, referenceConfigurationX[i]};
-    referenceConfiguration[i][_1] = Rotation<double,3>::identity();
+    referenceConfiguration[i][_1] = GFE::Rotation<double,3>::identity();
   }
 
   // Select the reference configuration as initial iterate
@@ -125,14 +125,14 @@ int main (int argc, char *argv[]) try
   FieldVector<double,3> axis = {1,0,0};
   double angle = 0;
 
-  x[rightBoundaryDof][_1] = Rotation<double,3>(axis, M_PI*angle/180);
+  x[rightBoundaryDof][_1] = GFE::Rotation<double,3>(axis, M_PI*angle/180);
 
   //////////////////////////////////////////////
   //  Create the energy and assembler
   //////////////////////////////////////////////
 
   using ATargetSpace = TargetSpace::rebind<adouble>::other;
-  using GeodesicInterpolationRule  = LocalGeodesicFEFunction<1, double, ScalarBasis::LocalView::Tree::FiniteElement, ATargetSpace>;
+  using GeodesicInterpolationRule  = GFE::LocalGeodesicFEFunction<1, double, ScalarBasis::LocalView::Tree::FiniteElement, ATargetSpace>;
   using ProjectedInterpolationRule = GFE::LocalProjectedFEFunction<1, double, ScalarBasis::LocalView::Tree::FiniteElement, ATargetSpace>;
 
   // Assembler using ADOL-C
@@ -157,16 +157,16 @@ int main (int argc, char *argv[]) try
   else
     DUNE_THROW(Exception, "Unknown interpolation method " << interpolationMethod << " requested!");
 
-  LocalGeodesicFEADOLCStiffness<ScalarBasis,
+  GFE::LocalGeodesicFEADOLCStiffness<ScalarBasis,
       TargetSpace> localStiffness(localRodEnergy);
 
-  GeodesicFEAssembler<ScalarBasis,TargetSpace> rodAssembler(gridView, localStiffness);
+  GFE::GeodesicFEAssembler<ScalarBasis,TargetSpace> rodAssembler(gridView, localStiffness);
 
   /////////////////////////////////////////////
   //   Create a solver for the rod problem
   /////////////////////////////////////////////
 
-  RiemannianTrustRegionSolver<ScalarBasis,TargetSpace> solver;
+  GFE::RiemannianTrustRegionSolver<ScalarBasis,TargetSpace> solver;
 
   solver.setup(grid,
                &rodAssembler,
diff --git a/test/embeddedglobalgfefunctiontest.cc b/test/embeddedglobalgfefunctiontest.cc
index 861b5846..2b798d5f 100644
--- a/test/embeddedglobalgfefunctiontest.cc
+++ b/test/embeddedglobalgfefunctiontest.cc
@@ -30,12 +30,12 @@ int main(int argc, char** argv)
   auto basis = Functions::BasisFactory::makeBasis(gridView, lagrange<2>());
 
   // Make a test coefficient set
-  using TargetSpace = UnitVector<double,3>;
+  using TargetSpace = GFE::UnitVector<double,3>;
 
   std::vector<TargetSpace> coefficients(basis.size());
   std::fill(coefficients.begin(), coefficients.end(), FieldVector<double,3>({1,0,0}));
 
-  using GeodesicInterpolationRule  = LocalGeodesicFEFunction<dim, double, decltype(basis)::LocalView::Tree::FiniteElement, TargetSpace>;
+  using GeodesicInterpolationRule = GFE::LocalGeodesicFEFunction<dim, double, decltype(basis)::LocalView::Tree::FiniteElement, TargetSpace>;
   GFE::EmbeddedGlobalGFEFunction<decltype(basis),GeodesicInterpolationRule,TargetSpace> testFunction(basis, coefficients);
 
   // Evaluate the function at the element centers
diff --git a/test/filmonsubstratetest.cc b/test/filmonsubstratetest.cc
index d3492855..5212bbff 100644
--- a/test/filmonsubstratetest.cc
+++ b/test/filmonsubstratetest.cc
@@ -143,9 +143,9 @@ int main (int argc, char *argv[])
 
   using namespace Dune::Indices;
 
-  typedef std::vector<RealTuple<double,dim> > DisplacementVector;
-  typedef std::vector<Rotation<double,dim> > RotationVector;
-  const int dimRotation = Rotation<double,dim>::embeddedDim;
+  typedef std::vector<GFE::RealTuple<double,dim> > DisplacementVector;
+  typedef std::vector<GFE::Rotation<double,dim> > RotationVector;
+  const int dimRotation = GFE::Rotation<double,dim>::embeddedDim;
   typedef TupleVector<DisplacementVector, RotationVector> SolutionType;
 
   /////////////////////////////////////////////////////////////
@@ -241,7 +241,7 @@ int main (int argc, char *argv[])
 #endif
 
   //Create BitVector matching the tangential space
-  const int dimRotationTangent = Rotation<double,dim>::TangentVector::dimension;
+  const int dimRotationTangent = GFE::Rotation<double,dim>::TangentVector::dimension;
   typedef MultiTypeBlockVector<std::vector<FieldVector<double,dim> >, std::vector<FieldVector<double,dimRotationTangent> > > VectorForBit;
   typedef Solvers::DefaultBitVector_t<VectorForBit> BitVector;
 
@@ -276,7 +276,7 @@ int main (int argc, char *argv[])
     x[_0][i] = identity[i];
 
   for (std::size_t i=0; i<x[_1].size(); ++i)
-    x[_1][i] = Rotation<double,dim>::identity();
+    x[_1][i] = GFE::Rotation<double,dim>::identity();
 
 
   /////////////////////////////////////////////////////////////
@@ -340,7 +340,7 @@ int main (int argc, char *argv[])
   materialParameters["b2"] = "1.0";
   materialParameters["b3"] = "1.0";
 
-  using ActiveRigidBodyMotion = GFE::ProductManifold<RealTuple<adouble,dim>, Rotation<adouble,dim> >;
+  using ActiveRigidBodyMotion = GFE::ProductManifold<GFE::RealTuple<adouble,dim>, GFE::Rotation<adouble,dim> >;
 
   std::shared_ptr<Elasticity::LocalDensity<dim,ValueType> > elasticDensity;
   elasticDensity = std::make_shared<Elasticity::MooneyRivlinDensity<dim,ValueType> >(materialParameters);
@@ -349,13 +349,13 @@ int main (int argc, char *argv[])
   auto elasticDensityWrapped = std::make_shared<GFE::DuneElasticityDensity<Element,ActiveRigidBodyMotion,0> >(elasticDensity);
 
   // Select which type of geometric interpolation to use
-  using LocalDeformationInterpolationRule = LocalGeodesicFEFunction<dim, GridType::ctype, decltype(deformationFEBasis.localView().tree().finiteElement()), RealTuple<adouble,dim> >;
-  using LocalOrientationInterpolationRule = LocalGeodesicFEFunction<dim, GridType::ctype, decltype(orientationFEBasis.localView().tree().finiteElement()), Rotation<adouble,dim> >;
+  using LocalDeformationInterpolationRule = GFE::LocalGeodesicFEFunction<dim, GridType::ctype, decltype(deformationFEBasis.localView().tree().finiteElement()), GFE::RealTuple<adouble,dim> >;
+  using LocalOrientationInterpolationRule = GFE::LocalGeodesicFEFunction<dim, GridType::ctype, decltype(orientationFEBasis.localView().tree().finiteElement()), GFE::Rotation<adouble,dim> >;
 
   using LocalInterpolationRule = std::tuple<LocalDeformationInterpolationRule,LocalOrientationInterpolationRule>;
 
   auto elasticEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, LocalInterpolationRule, ActiveRigidBodyMotion> >(elasticDensityWrapped);
-  auto neumannEnergy = std::make_shared<GFE::NeumannEnergy<CompositeBasis, RealTuple<ValueType,targetDim>, Rotation<ValueType,dim> > >(neumannBoundary,neumannFunction);
+  auto neumannEnergy = std::make_shared<GFE::NeumannEnergy<CompositeBasis, GFE::RealTuple<ValueType,targetDim>, GFE::Rotation<ValueType,dim> > >(neumannBoundary,neumannFunction);
 
   using Intersection = typename GridView::Intersection;
   auto cosseratShellDensity = std::make_shared<GFE::CosseratShellDensity<Intersection, adouble> >(
@@ -369,28 +369,28 @@ int main (int argc, char *argv[])
     &surfaceShellBoundary,
     stressFreeShellFunction);
 
-  using RBM = GFE::ProductManifold<RealTuple<double, dim>,Rotation<double,dim> >;
+  using RBM = GFE::ProductManifold<GFE::RealTuple<double, dim>,GFE::Rotation<double,dim> >;
 
-  auto sumEnergy = std::make_shared<GFE::SumEnergy<CompositeBasis, RealTuple<ValueType,targetDim>, Rotation<ValueType,targetDim> > >();
+  auto sumEnergy = std::make_shared<GFE::SumEnergy<CompositeBasis, GFE::RealTuple<ValueType,targetDim>, GFE::Rotation<ValueType,targetDim> > >();
   sumEnergy->addLocalEnergy(neumannEnergy);
   sumEnergy->addLocalEnergy(elasticEnergy);
   sumEnergy->addLocalEnergy(surfaceCosseratEnergy);
 
-  LocalGeodesicFEADOLCStiffness<CompositeBasis,RBM> localGFEADOLCStiffness(sumEnergy);
-  MixedGFEAssembler<CompositeBasis,RBM> mixedAssembler(compositeBasis, localGFEADOLCStiffness);
+  GFE::LocalGeodesicFEADOLCStiffness<CompositeBasis,RBM> localGFEADOLCStiffness(sumEnergy);
+  GFE::MixedGFEAssembler<CompositeBasis,RBM> mixedAssembler(compositeBasis, localGFEADOLCStiffness);
 
   ////////////////////////////////////////////////////////
   //   Set Dirichlet values
   ////////////////////////////////////////////////////////
 
-  BitSetVector<RealTuple<double,dim>::TangentVector::dimension> deformationDirichletDofs(dirichletDofs[_0].size(), false);
+  BitSetVector<GFE::RealTuple<double,dim>::TangentVector::dimension> deformationDirichletDofs(dirichletDofs[_0].size(), false);
   for (std::size_t i = 0; i < dirichletDofs[_0].size(); i++)
-    for (int j = 0; j < RealTuple<double,dim>::TangentVector::dimension; j++)
+    for (int j = 0; j < GFE::RealTuple<double,dim>::TangentVector::dimension; j++)
       deformationDirichletDofs[i][j] = dirichletDofs[_0][i][j];
 
-  BitSetVector<Rotation<double,dim>::TangentVector::dimension> orientationDirichletDofs(dirichletDofs[_1].size(), false);
+  BitSetVector<GFE::Rotation<double,dim>::TangentVector::dimension> orientationDirichletDofs(dirichletDofs[_1].size(), false);
   for (std::size_t i = 0; i < dirichletDofs[_1].size(); i++)
-    for (int j = 0; j < Rotation<double,dim>::TangentVector::dimension; j++)
+    for (int j = 0; j < GFE::Rotation<double,dim>::TangentVector::dimension; j++)
       orientationDirichletDofs[i][j] = dirichletDofs[_1][i][j];
 
 #if !MIXED_SPACE
@@ -421,7 +421,7 @@ int main (int argc, char *argv[])
   if (solverType == "trustRegion")
   {
 #if MIXED_SPACE
-    MixedRiemannianTrustRegionSolver<GridType, CompositeBasis, DeformationFEBasis, RealTuple<double,dim>, OrientationFEBasis, Rotation<double,dim> > solver;
+    GFE::MixedRiemannianTrustRegionSolver<GridType, CompositeBasis, DeformationFEBasis, GFE::RealTuple<double,dim>, OrientationFEBasis, GFE::Rotation<double,dim> > solver;
     solver.setup(*grid,
                  &mixedAssembler,
                  deformationFEBasis,
@@ -444,7 +444,7 @@ int main (int argc, char *argv[])
     solver.solve();
     x = solver.getSol();
 #else
-    RiemannianTrustRegionSolver<DeformationFEBasis, RBM, GFEAssemblerWrapper> solver;
+    GFE::RiemannianTrustRegionSolver<DeformationFEBasis, RBM, GFEAssemblerWrapper> solver;
     solver.setup(*grid,
                  &assembler,
                  xRBM,
@@ -473,7 +473,7 @@ int main (int argc, char *argv[])
   } else {    // solverType == "proximalNewton"
 
 #if MIXED_SPACE
-    GFE::MixedRiemannianProximalNewtonSolver<CompositeBasis, DeformationFEBasis, RealTuple<double,dim>, OrientationFEBasis, Rotation<double,dim>, BitVector> solver;
+    GFE::MixedRiemannianProximalNewtonSolver<CompositeBasis, DeformationFEBasis, GFE::RealTuple<double,dim>, OrientationFEBasis, GFE::Rotation<double,dim>, BitVector> solver;
     solver.setup(*grid,
                  &mixedAssembler,
                  x,
@@ -486,7 +486,7 @@ int main (int argc, char *argv[])
     solver.solve();
     x = solver.getSol();
 #else
-    RiemannianProximalNewtonSolver<DeformationFEBasis, RBM, GFEAssemblerWrapper> solver;
+    GFE::RiemannianProximalNewtonSolver<DeformationFEBasis, RBM, GFEAssemblerWrapper> solver;
     solver.setup(*grid,
                  &assembler,
                  xRBM,
diff --git a/test/harmonicmaptest.cc b/test/harmonicmaptest.cc
index 580051e1..5f1c0446 100644
--- a/test/harmonicmaptest.cc
+++ b/test/harmonicmaptest.cc
@@ -35,7 +35,7 @@
 const int dim = 2;
 
 // Image space of the geodesic fe functions
-typedef UnitVector<double,3> TargetSpace;
+typedef Dune::GFE::UnitVector<double,3> TargetSpace;
 
 const int order = 1;
 
@@ -150,7 +150,7 @@ int main (int argc, char *argv[])
 
 
 #if GEODESICINTERPOLATION
-  using InterpolationRule = LocalGeodesicFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, TargetSpace>;
+  using InterpolationRule = GFE::LocalGeodesicFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, TargetSpace>;
 #else
 #if CONFORMING
   using InterpolationRule = GFE::LocalProjectedFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, TargetSpace,true>;
@@ -163,13 +163,13 @@ int main (int argc, char *argv[])
 
   GFE::LocalIntegralStiffness<FEBasis,InterpolationRule,TargetSpace> localGFEADOLCStiffness(harmonicDensity);
 
-  GeodesicFEAssembler<FEBasis,TargetSpace> assembler(feBasis, localGFEADOLCStiffness);
+  GFE::GeodesicFEAssembler<FEBasis,TargetSpace> assembler(feBasis, localGFEADOLCStiffness);
 
   ///////////////////////////////////////////////////
   //  Create a Riemannian trust-region solver
   ///////////////////////////////////////////////////
 
-  RiemannianTrustRegionSolver<FEBasis,TargetSpace> solver;
+  GFE::RiemannianTrustRegionSolver<FEBasis,TargetSpace> solver;
   solver.setup(*grid,
                &assembler,
                x,
diff --git a/test/interpolationderivativestest.cc b/test/interpolationderivativestest.cc
index f64a5074..d633f78e 100644
--- a/test/interpolationderivativestest.cc
+++ b/test/interpolationderivativestest.cc
@@ -391,7 +391,7 @@ TestSuite checkDerivatives()
   Functions::LagrangeBasis<GridView,order> scalarBasis(gridView);
 
   std::vector<TargetSpace> testPoints;
-  ValueFactory<TargetSpace>::get(testPoints);
+  GFE::ValueFactory<TargetSpace>::get(testPoints);
 
   // TODO: Make sure the list of test points is longer than this.
   const std::size_t nDofs = scalarBasis.dimension();
@@ -405,7 +405,7 @@ TestSuite checkDerivatives()
   /////////////////////////////////////////////////////////////////////////
 
   // Define the two possible interpolation rules
-  using GeodesicInterpolationRule = LocalGeodesicFEFunction<domainDim,
+  using GeodesicInterpolationRule = GFE::LocalGeodesicFEFunction<domainDim,
       typename Grid::ctype,
       decltype(scalarBasis.localView().tree().finiteElement()),
       TargetSpace>;
@@ -556,16 +556,16 @@ int main (int argc, char *argv[])
 
   // Test the UnitSphere class and geodesic interpolation.
   // This uses the default derivatives implementation (using ADOL-C)
-  test.subTest(checkDerivatives<UnitVector<double,3>, InterpolationType::Geodesic >());
+  test.subTest(checkDerivatives<GFE::UnitVector<double,3>, InterpolationType::Geodesic >());
 
   // Test the RealTuple class, both with geodesic and projection-based interpolation
   // Both are specialized
-  test.subTest(checkDerivatives<RealTuple<double,3>, InterpolationType::Geodesic>());
-  test.subTest(checkDerivatives<RealTuple<double,3>, InterpolationType::ProjectionBased>());
+  test.subTest(checkDerivatives<GFE::RealTuple<double,3>, InterpolationType::Geodesic>());
+  test.subTest(checkDerivatives<GFE::RealTuple<double,3>, InterpolationType::ProjectionBased>());
 
   // Test the UnitVector class with projection-based interpolation
   // This is also specialized.
-  test.subTest(checkDerivatives<UnitVector<double,3>, InterpolationType::ProjectionBased>());
+  test.subTest(checkDerivatives<GFE::UnitVector<double,3>, InterpolationType::ProjectionBased>());
 
   return test.exit();
 }
diff --git a/test/localadolcstiffnesstest.cc b/test/localadolcstiffnesstest.cc
index ee545ad2..946cbf86 100644
--- a/test/localadolcstiffnesstest.cc
+++ b/test/localadolcstiffnesstest.cc
@@ -52,7 +52,7 @@ using namespace Dune;
 const int dim = 2;
 
 // Image space of the geodesic fe functions
-using TargetSpace = GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> >;
+using TargetSpace = GFE::ProductManifold<GFE::RealTuple<double,3>,GFE::Rotation<double,3> >;
 
 // Compare two matrices
 void compareMatrices(const Matrix<double>& matrixA, std::string nameA,
@@ -112,7 +112,7 @@ int main (int argc, char *argv[]) try
 
   using namespace Functions::BasisFactory;
 
-  const int dimRotation = Rotation<double,3>::TangentVector::dimension;
+  const int dimRotation = GFE::Rotation<double,3>::TangentVector::dimension;
 
   auto tangentBasis = makeBasis(
     gridView,
@@ -187,14 +187,14 @@ int main (int argc, char *argv[]) try
   using ATargetSpace = typename TargetSpace::template rebind<adouble>::other;
 
   // Select geometric finite element interpolation method
-  using AInterpolationRule = LocalGeodesicFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, ATargetSpace>;
+  using AInterpolationRule = GFE::LocalGeodesicFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, ATargetSpace>;
 
   auto activeDensity = std::make_shared<GFE::PlanarCosseratShellDensity<GridType::Codim<0>::Entity, adouble> >(materialParameters);
 
   auto activeCosseratLocalEnergy = std::make_shared<GFE::LocalIntegralEnergy<TangentBasis,AInterpolationRule,ATargetSpace> >(activeDensity);
 
   // The actual assembler
-  LocalGeodesicFEADOLCStiffness<TangentBasis,
+  GFE::LocalGeodesicFEADOLCStiffness<TangentBasis,
       TargetSpace> localGFEADOLCStiffness(activeCosseratLocalEnergy);
 
   //////////////////////////////////////////////////////
@@ -202,14 +202,14 @@ int main (int argc, char *argv[]) try
   //////////////////////////////////////////////////////
 
   // Select geometric finite element interpolation method
-  using InterpolationRule = LocalGeodesicFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, TargetSpace>;
+  using InterpolationRule = GFE::LocalGeodesicFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, TargetSpace>;
 
   auto cosseratDensity = std::make_shared<GFE::PlanarCosseratShellDensity<GridType::Codim<0>::Entity, double> >(materialParameters);
 
   auto cosseratLocalEnergy = std::make_shared<GFE::LocalIntegralEnergy<TangentBasis,InterpolationRule,TargetSpace> >(cosseratDensity);
 
   // The actual assembler
-  LocalGeodesicFEFDStiffness<TangentBasis,
+  GFE::LocalGeodesicFEFDStiffness<TangentBasis,
       TargetSpace,
       FDType> localGFEFDStiffness(cosseratLocalEnergy.get());
 
diff --git a/test/localgeodesicfefunctiontest.cc b/test/localgeodesicfefunctiontest.cc
index 1f02784f..db8c1544 100644
--- a/test/localgeodesicfefunctiontest.cc
+++ b/test/localgeodesicfefunctiontest.cc
@@ -87,9 +87,9 @@ void testPermutationInvariance(const std::vector<TargetSpace>& corners)
   cornersRotated1[1] = cornersRotated2[0] = corners[2];
   cornersRotated1[2] = cornersRotated2[1] = corners[0];
 
-  LocalGeodesicFEFunction<2,double,LocalFiniteElement,TargetSpace> f0(feCache.get(simplex), corners);
-  LocalGeodesicFEFunction<2,double,LocalFiniteElement,TargetSpace> f1(feCache.get(simplex), cornersRotated1);
-  LocalGeodesicFEFunction<2,double,LocalFiniteElement,TargetSpace> f2(feCache.get(simplex), cornersRotated2);
+  GFE::LocalGeodesicFEFunction<2,double,LocalFiniteElement,TargetSpace> f0(feCache.get(simplex), corners);
+  GFE::LocalGeodesicFEFunction<2,double,LocalFiniteElement,TargetSpace> f1(feCache.get(simplex), cornersRotated1);
+  GFE::LocalGeodesicFEFunction<2,double,LocalFiniteElement,TargetSpace> f2(feCache.get(simplex), cornersRotated2);
 
   // A quadrature rule as a set of test points
   int quadOrder = 3;
@@ -124,7 +124,7 @@ void testPermutationInvariance(const std::vector<TargetSpace>& corners)
 }
 
 template <int domainDim, class TargetSpace>
-void testDerivative(const LocalGeodesicFEFunction<domainDim,double,typename LagrangeLocalFiniteElementCache<double,double,domainDim,1>::FiniteElementType, TargetSpace>& f)
+void testDerivative(const GFE::LocalGeodesicFEFunction<domainDim,double,typename LagrangeLocalFiniteElementCache<double,double,domainDim,1>::FiniteElementType, TargetSpace>& f)
 {
   static const int embeddedDim = TargetSpace::EmbeddedTangentVector::dimension;
 
@@ -158,7 +158,7 @@ void testDerivative(const LocalGeodesicFEFunction<domainDim,double,typename Lagr
 
 
 template <int domainDim, class TargetSpace>
-void testDerivativeOfValueWRTCoefficients(const LocalGeodesicFEFunction<domainDim,double,typename LagrangeLocalFiniteElementCache<double,double,domainDim,1>::FiniteElementType, TargetSpace>& f)
+void testDerivativeOfValueWRTCoefficients(const GFE::LocalGeodesicFEFunction<domainDim,double,typename LagrangeLocalFiniteElementCache<double,double,domainDim,1>::FiniteElementType, TargetSpace>& f)
 {
   static const int embeddedDim = TargetSpace::EmbeddedTangentVector::dimension;
 
@@ -202,7 +202,7 @@ void testDerivativeOfValueWRTCoefficients(const LocalGeodesicFEFunction<domainDi
 }
 
 template <int domainDim, class TargetSpace>
-void testDerivativeOfGradientWRTCoefficients(const LocalGeodesicFEFunction<domainDim,double,typename LagrangeLocalFiniteElementCache<double,double,domainDim,1>::FiniteElementType, TargetSpace>& f)
+void testDerivativeOfGradientWRTCoefficients(const GFE::LocalGeodesicFEFunction<domainDim,double,typename LagrangeLocalFiniteElementCache<double,double,domainDim,1>::FiniteElementType, TargetSpace>& f)
 {
   static const int embeddedDim = TargetSpace::EmbeddedTangentVector::dimension;
 
@@ -220,11 +220,11 @@ void testDerivativeOfGradientWRTCoefficients(const LocalGeodesicFEFunction<domai
     for (size_t i=0; i<f.size(); i++) {
 
       // evaluate actual derivative
-      Tensor3<double, embeddedDim, embeddedDim, domainDim> derivative;
+      GFE::Tensor3<double, embeddedDim, embeddedDim, domainDim> derivative;
       f.evaluateDerivativeOfGradientWRTCoefficient(quadPos, i, derivative);
 
       // evaluate fd approximation of derivative
-      Tensor3<double, embeddedDim, embeddedDim, domainDim> fdDerivative;
+      GFE::Tensor3<double, embeddedDim, embeddedDim, domainDim> fdDerivative;
       f.evaluateFDDerivativeOfGradientWRTCoefficient(quadPos, i, fdDerivative);
 
       if ( (derivative - fdDerivative).infinity_norm() > eps ) {
@@ -252,7 +252,7 @@ void test(const GeometryType& element)
   std::cout << " --- Testing " << className<TargetSpace>() << ", domain dimension: " << element.dim() << " ---" << std::endl;
 
   std::vector<TargetSpace> testPoints;
-  ValueFactory<TargetSpace>::get(testPoints);
+  GFE::ValueFactory<TargetSpace>::get(testPoints);
 
   int nTestPoints = testPoints.size();
   size_t nVertices = Dune::ReferenceElements<double,domainDim>::general(element).size(domainDim);
@@ -260,7 +260,7 @@ void test(const GeometryType& element)
   // Set up elements of the target space
   std::vector<TargetSpace> corners(nVertices);
 
-  MultiIndex index(nVertices, nTestPoints);
+  GFE::MultiIndex index(nVertices, nTestPoints);
   int numIndices = index.cycle();
 
   for (int i=0; i<numIndices; i++, ++index) {
@@ -275,7 +275,7 @@ void test(const GeometryType& element)
     LagrangeLocalFiniteElementCache<double,double,domainDim,1> feCache;
     typedef typename LagrangeLocalFiniteElementCache<double,double,domainDim,1>::FiniteElementType LocalFiniteElement;
 
-    LocalGeodesicFEFunction<domainDim,double,LocalFiniteElement,TargetSpace> f(feCache.get(element),corners);
+    GFE::LocalGeodesicFEFunction<domainDim,double,LocalFiniteElement,TargetSpace> f(feCache.get(element),corners);
 
     //testPermutationInvariance(corners);
     testDerivative<domainDim>(f);
@@ -299,33 +299,33 @@ int main()
   //  Test functions on 1d elements
   ////////////////////////////////////////////////////////////////
 
-  test<RealTuple<double,1>,1>(GeometryTypes::simplex(1));
-  test<UnitVector<double,2>,1>(GeometryTypes::simplex(1));
-  test<UnitVector<double,3>,1>(GeometryTypes::simplex(1));
-  test<Rotation<double,3>,1>(GeometryTypes::simplex(1));
-  typedef Dune::GFE::ProductManifold<RealTuple<double,1>,Rotation<double,3>,UnitVector<double,2> > CrazyManifold;
+  test<GFE::RealTuple<double,1>,1>(GeometryTypes::simplex(1));
+  test<GFE::UnitVector<double,2>,1>(GeometryTypes::simplex(1));
+  test<GFE::UnitVector<double,3>,1>(GeometryTypes::simplex(1));
+  test<GFE::Rotation<double,3>,1>(GeometryTypes::simplex(1));
+  typedef GFE::ProductManifold<GFE::RealTuple<double,1>,GFE::Rotation<double,3>,GFE::UnitVector<double,2> > CrazyManifold;
   test<CrazyManifold,1>(GeometryTypes::simplex(1));
 
   ////////////////////////////////////////////////////////////////
   //  Test functions on 2d simplex elements
   ////////////////////////////////////////////////////////////////
 
-  test<RealTuple<double,1>,2>(GeometryTypes::simplex(2));
-  test<UnitVector<double,2>,2>(GeometryTypes::simplex(2));
-  test<UnitVector<double,3>,2>(GeometryTypes::simplex(2));
-  test<Rotation<double,3>,2>(GeometryTypes::simplex(2));
-  typedef Dune::GFE::ProductManifold<RealTuple<double,1>,Rotation<double,3>,UnitVector<double,2> > CrazyManifold;
+  test<GFE::RealTuple<double,1>,2>(GeometryTypes::simplex(2));
+  test<GFE::UnitVector<double,2>,2>(GeometryTypes::simplex(2));
+  test<GFE::UnitVector<double,3>,2>(GeometryTypes::simplex(2));
+  test<GFE::Rotation<double,3>,2>(GeometryTypes::simplex(2));
+  typedef GFE::ProductManifold<GFE::RealTuple<double,1>,GFE::Rotation<double,3>,GFE::UnitVector<double,2> > CrazyManifold;
   test<CrazyManifold,2>(GeometryTypes::simplex(2));
 
   ////////////////////////////////////////////////////////////////
   //  Test functions on 2d quadrilateral elements
   ////////////////////////////////////////////////////////////////
 
-  test<RealTuple<double,1>,2>(GeometryTypes::cube(2));
-  test<UnitVector<double,2>,2>(GeometryTypes::cube(2));
-  test<UnitVector<double,3>,2>(GeometryTypes::cube(2));
-  test<Rotation<double,3>,2>(GeometryTypes::cube(2));
-  typedef Dune::GFE::ProductManifold<RealTuple<double,1>,Rotation<double,3>,UnitVector<double,2> > CrazyManifold;
+  test<GFE::RealTuple<double,1>,2>(GeometryTypes::cube(2));
+  test<GFE::UnitVector<double,2>,2>(GeometryTypes::cube(2));
+  test<GFE::UnitVector<double,3>,2>(GeometryTypes::cube(2));
+  test<GFE::Rotation<double,3>,2>(GeometryTypes::cube(2));
+  typedef GFE::ProductManifold<GFE::RealTuple<double,1>,GFE::Rotation<double,3>,GFE::UnitVector<double,2> > CrazyManifold;
   test<CrazyManifold,2>(GeometryTypes::cube(2));
 
 }
diff --git a/test/localintegralstiffnesstest.cc b/test/localintegralstiffnesstest.cc
index 55691ac5..37b4b6f9 100644
--- a/test/localintegralstiffnesstest.cc
+++ b/test/localintegralstiffnesstest.cc
@@ -45,7 +45,7 @@ enum InterpolationType {Geodesic, ProjectionBased, Nonconforming};
 template <class GridView, InterpolationType interpolationType>
 int testHarmonicMapIntoSphere(TestSuite& test, const GridView& gridView)
 {
-  using TargetSpace = UnitVector<double,dim>;
+  using TargetSpace = GFE::UnitVector<double,dim>;
 
   // Finite element order
   const int order = 1;
@@ -97,8 +97,8 @@ int testHarmonicMapIntoSphere(TestSuite& test, const GridView& gridView)
   //  Set up the two assemblers
   //////////////////////////////////////////////////////////////
 
-  std::shared_ptr<GeodesicFEAssembler<FEBasis,TargetSpace> > assemblerADOLC;
-  std::shared_ptr<LocalGeodesicFEStiffness<FEBasis,TargetSpace> > localIntegralStiffness;
+  std::shared_ptr<GFE::GeodesicFEAssembler<FEBasis,TargetSpace> > assemblerADOLC;
+  std::shared_ptr<GFE::LocalGeodesicFEStiffness<FEBasis,TargetSpace> > localIntegralStiffness;
 
   using Element = typename GridView::template Codim<0>::Entity;
   auto harmonicDensity = std::make_shared<GFE::HarmonicDensity<Element,TargetSpace> >();
@@ -107,12 +107,12 @@ int testHarmonicMapIntoSphere(TestSuite& test, const GridView& gridView)
   if constexpr (interpolationType==Geodesic)
   {
     std::cout << "Using geodesic interpolation" << std::endl;
-    using LocalInterpolationRule = LocalGeodesicFEFunction<dim,
+    using LocalInterpolationRule = GFE::LocalGeodesicFEFunction<dim,
         typename GridView::ctype,
         decltype(feBasis.localView().tree().finiteElement()),
         TargetSpace>;
 
-    using ALocalInterpolationRule = LocalGeodesicFEFunction<dim,
+    using ALocalInterpolationRule = GFE::LocalGeodesicFEFunction<dim,
         typename GridView::ctype,
         decltype(feBasis.localView().tree().finiteElement()),
         ATargetSpace>;
@@ -120,8 +120,8 @@ int testHarmonicMapIntoSphere(TestSuite& test, const GridView& gridView)
     // Assemble using the old assembler
     auto energy = std::make_shared<GFE::LocalIntegralEnergy<FEBasis, ALocalInterpolationRule,ATargetSpace> >(harmonicDensityA);
 
-    auto localGFEADOLCStiffness = std::make_shared<LocalGeodesicFEADOLCStiffness<FEBasis,TargetSpace> >(energy);
-    assemblerADOLC = std::make_shared<GeodesicFEAssembler<FEBasis,TargetSpace> >(feBasis, localGFEADOLCStiffness);
+    auto localGFEADOLCStiffness = std::make_shared<GFE::LocalGeodesicFEADOLCStiffness<FEBasis,TargetSpace> >(energy);
+    assemblerADOLC = std::make_shared<GFE::GeodesicFEAssembler<FEBasis,TargetSpace> >(feBasis, localGFEADOLCStiffness);
 
     // Assemble using the new assembler
     localIntegralStiffness = std::make_shared<GFE::LocalIntegralStiffness<FEBasis,LocalInterpolationRule,TargetSpace> >(harmonicDensity);
@@ -146,15 +146,15 @@ int testHarmonicMapIntoSphere(TestSuite& test, const GridView& gridView)
     // Assemble using the old assembler
     auto energy = std::make_shared<GFE::LocalIntegralEnergy<FEBasis, ALocalInterpolationRule,ATargetSpace> >(harmonicDensityA);
 
-    auto localGFEADOLCStiffness = std::make_shared<LocalGeodesicFEADOLCStiffness<FEBasis,TargetSpace> >(energy);
-    assemblerADOLC = std::make_shared<GeodesicFEAssembler<FEBasis,TargetSpace> >(feBasis, localGFEADOLCStiffness);
+    auto localGFEADOLCStiffness = std::make_shared<GFE::LocalGeodesicFEADOLCStiffness<FEBasis,TargetSpace> >(energy);
+    assemblerADOLC = std::make_shared<GFE::GeodesicFEAssembler<FEBasis,TargetSpace> >(feBasis, localGFEADOLCStiffness);
 
     // Assemble using the new assembler
     localIntegralStiffness = std::make_shared<GFE::LocalIntegralStiffness<FEBasis,LocalInterpolationRule,TargetSpace> >(harmonicDensity);
   }
 
   // TODO: The assembler should really get the tangent basis.
-  auto assemblerIntegralStiffness = std::make_shared<GeodesicFEAssembler<FEBasis,TargetSpace> >(feBasis, localIntegralStiffness);
+  auto assemblerIntegralStiffness = std::make_shared<GFE::GeodesicFEAssembler<FEBasis,TargetSpace> >(feBasis, localIntegralStiffness);
 
   //////////////////////////////////////////////////////////////
   //  Assemble
@@ -224,8 +224,8 @@ int testCosseratBulkModel(TestSuite& test, const GridView& gridView)
   const int deformationOrder = 2;
   const int rotationOrder = 1;
 
-  const static int deformationBlocksize = RealTuple<double,dim>::TangentVector::dimension;
-  const static int orientationBlocksize = Rotation<double,dim>::TangentVector::dimension;
+  const static int deformationBlocksize = GFE::RealTuple<double,dim>::TangentVector::dimension;
+  const static int orientationBlocksize = GFE::Rotation<double,dim>::TangentVector::dimension;
 
   using CorrectionType0 = BlockVector<FieldVector<double, deformationBlocksize> >;
   using CorrectionType1 = BlockVector<FieldVector<double, orientationBlocksize> >;
@@ -242,7 +242,7 @@ int testCosseratBulkModel(TestSuite& test, const GridView& gridView)
   //  Construct all needed function space bases
   //////////////////////////////////////////////////////////////////////////////////
 
-  const int dimRotation = Rotation<double,dim>::embeddedDim;
+  const int dimRotation = GFE::Rotation<double,dim>::embeddedDim;
   auto compositeBasis = makeBasis(
     gridView,
     composite(
@@ -260,7 +260,7 @@ int testCosseratBulkModel(TestSuite& test, const GridView& gridView)
   //  Construct the configuration where to assemble the tangent problem
   /////////////////////////////////////////////////////////////////////////
 
-  using CoefficientVector = TupleVector<std::vector<RealTuple<double,dim> >,std::vector<Rotation<double,dim> > >;
+  using CoefficientVector = TupleVector<std::vector<GFE::RealTuple<double,dim> >,std::vector<GFE::Rotation<double,dim> > >;
   CoefficientVector x;
   x[_0].resize(compositeBasis.size({0}));
   x[_1].resize(compositeBasis.size({1}));
@@ -279,7 +279,7 @@ int testCosseratBulkModel(TestSuite& test, const GridView& gridView)
     return y;
   });
   for (size_t i = 0; i < x[_1].size(); i++)
-    x[_1][i] = Rotation<double,dim>(initialIterate[_1][i]);
+    x[_1][i] = GFE::Rotation<double,dim>(initialIterate[_1][i]);
 
   // Set of material parameters
   ParameterTree parameters;
@@ -295,8 +295,8 @@ int testCosseratBulkModel(TestSuite& test, const GridView& gridView)
   parameters["b2"] = "1";
   parameters["b3"] = "1";
 
-  using RigidBodyMotion  = GFE::ProductManifold<RealTuple<double,dim>, Rotation<double,dim> >;
-  using ARigidBodyMotion = GFE::ProductManifold<RealTuple<adouble,dim>, Rotation<adouble,dim> >;
+  using RigidBodyMotion  = GFE::ProductManifold<GFE::RealTuple<double,dim>, GFE::Rotation<double,dim> >;
+  using ARigidBodyMotion = GFE::ProductManifold<GFE::RealTuple<adouble,dim>, GFE::Rotation<adouble,dim> >;
 
   using Element = typename GridView::template Codim<0>::Entity;
 
@@ -314,29 +314,29 @@ int testCosseratBulkModel(TestSuite& test, const GridView& gridView)
   //  Set up the two assemblers
   //////////////////////////////////////////////////////////////
 
-  std::shared_ptr<MixedGFEAssembler<CompositeBasis,RigidBodyMotion> > assemblerADOLC;
-  std::shared_ptr<LocalGeodesicFEStiffness<CompositeBasis,RigidBodyMotion> > localIntegralStiffness;
+  std::shared_ptr<GFE::MixedGFEAssembler<CompositeBasis,RigidBodyMotion> > assemblerADOLC;
+  std::shared_ptr<GFE::LocalGeodesicFEStiffness<CompositeBasis,RigidBodyMotion> > localIntegralStiffness;
 
   if constexpr (interpolationType==Geodesic)
   {
     std::cout << "Using geodesic interpolation" << std::endl;
 
-    using LocalDeformationInterpolationRule = LocalGeodesicFEFunction<dim, typename GridView::ctype, decltype(deformationFEBasis.localView().tree().finiteElement()), RealTuple<double,dim> >;
-    using LocalOrientationInterpolationRule = LocalGeodesicFEFunction<dim, typename GridView::ctype, decltype(orientationFEBasis.localView().tree().finiteElement()), Rotation<double,dim> >;
+    using LocalDeformationInterpolationRule = GFE::LocalGeodesicFEFunction<dim, typename GridView::ctype, decltype(deformationFEBasis.localView().tree().finiteElement()), GFE::RealTuple<double,dim> >;
+    using LocalOrientationInterpolationRule = GFE::LocalGeodesicFEFunction<dim, typename GridView::ctype, decltype(orientationFEBasis.localView().tree().finiteElement()), GFE::Rotation<double,dim> >;
 
     using LocalInterpolationRule = std::tuple<LocalDeformationInterpolationRule,LocalOrientationInterpolationRule>;
 
-    using ALocalDeformationInterpolationRule = LocalGeodesicFEFunction<dim, typename GridView::ctype, decltype(deformationFEBasis.localView().tree().finiteElement()), RealTuple<adouble,dim> >;
-    using ALocalOrientationInterpolationRule = LocalGeodesicFEFunction<dim, typename GridView::ctype, decltype(orientationFEBasis.localView().tree().finiteElement()), Rotation<adouble,dim> >;
+    using ALocalDeformationInterpolationRule = GFE::LocalGeodesicFEFunction<dim, typename GridView::ctype, decltype(deformationFEBasis.localView().tree().finiteElement()), GFE::RealTuple<adouble,dim> >;
+    using ALocalOrientationInterpolationRule = GFE::LocalGeodesicFEFunction<dim, typename GridView::ctype, decltype(orientationFEBasis.localView().tree().finiteElement()), GFE::Rotation<adouble,dim> >;
 
     using ALocalInterpolationRule = std::tuple<ALocalDeformationInterpolationRule,ALocalOrientationInterpolationRule>;
 
     // Assemble using the ADOL-C assembler
     auto energy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, ALocalInterpolationRule,ARigidBodyMotion> >(aBulkCosseratDensity);
 
-    auto localGFEADOLCStiffness = std::make_shared<LocalGeodesicFEADOLCStiffness<CompositeBasis,RigidBodyMotion> >(energy);
+    auto localGFEADOLCStiffness = std::make_shared<GFE::LocalGeodesicFEADOLCStiffness<CompositeBasis,RigidBodyMotion> >(energy);
 
-    assemblerADOLC = std::make_shared<MixedGFEAssembler<CompositeBasis,RigidBodyMotion> >(compositeBasis, localGFEADOLCStiffness);
+    assemblerADOLC = std::make_shared<GFE::MixedGFEAssembler<CompositeBasis,RigidBodyMotion> >(compositeBasis, localGFEADOLCStiffness);
 
     // Assemble using the new assembler
     localIntegralStiffness = std::make_shared<GFE::LocalIntegralStiffness<CompositeBasis,LocalInterpolationRule,RigidBodyMotion> >(bulkCosseratDensity);
@@ -345,28 +345,28 @@ int testCosseratBulkModel(TestSuite& test, const GridView& gridView)
   {
     std::cout << "Using projection-based interpolation" << std::endl;
 
-    using LocalDeformationInterpolationRule = GFE::LocalProjectedFEFunction<dim, typename GridView::ctype, decltype(deformationFEBasis.localView().tree().finiteElement()), RealTuple<double,dim> >;
-    using LocalOrientationInterpolationRule = GFE::LocalProjectedFEFunction<dim, typename GridView::ctype, decltype(orientationFEBasis.localView().tree().finiteElement()), Rotation<double,dim> >;
+    using LocalDeformationInterpolationRule = GFE::LocalProjectedFEFunction<dim, typename GridView::ctype, decltype(deformationFEBasis.localView().tree().finiteElement()), GFE::RealTuple<double,dim> >;
+    using LocalOrientationInterpolationRule = GFE::LocalProjectedFEFunction<dim, typename GridView::ctype, decltype(orientationFEBasis.localView().tree().finiteElement()), GFE::Rotation<double,dim> >;
 
     using LocalInterpolationRule = std::tuple<LocalDeformationInterpolationRule,LocalOrientationInterpolationRule>;
 
-    using ALocalDeformationInterpolationRule = GFE::LocalProjectedFEFunction<dim, typename GridView::ctype, decltype(deformationFEBasis.localView().tree().finiteElement()), RealTuple<adouble,dim> >;
-    using ALocalOrientationInterpolationRule = GFE::LocalProjectedFEFunction<dim, typename GridView::ctype, decltype(orientationFEBasis.localView().tree().finiteElement()), Rotation<adouble,dim> >;
+    using ALocalDeformationInterpolationRule = GFE::LocalProjectedFEFunction<dim, typename GridView::ctype, decltype(deformationFEBasis.localView().tree().finiteElement()), GFE::RealTuple<adouble,dim> >;
+    using ALocalOrientationInterpolationRule = GFE::LocalProjectedFEFunction<dim, typename GridView::ctype, decltype(orientationFEBasis.localView().tree().finiteElement()), GFE::Rotation<adouble,dim> >;
 
     using ALocalInterpolationRule = std::tuple<ALocalDeformationInterpolationRule,ALocalOrientationInterpolationRule>;
 
     // Assemble using the ADOL-C assembler
     auto energy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, ALocalInterpolationRule,ARigidBodyMotion> >(aBulkCosseratDensity);
 
-    auto localGFEADOLCStiffness = std::make_shared<LocalGeodesicFEADOLCStiffness<CompositeBasis,RigidBodyMotion> >(energy);
+    auto localGFEADOLCStiffness = std::make_shared<GFE::LocalGeodesicFEADOLCStiffness<CompositeBasis,RigidBodyMotion> >(energy);
 
-    assemblerADOLC = std::make_shared<MixedGFEAssembler<CompositeBasis,RigidBodyMotion> >(compositeBasis, localGFEADOLCStiffness);
+    assemblerADOLC = std::make_shared<GFE::MixedGFEAssembler<CompositeBasis,RigidBodyMotion> >(compositeBasis, localGFEADOLCStiffness);
 
     // Assemble using the new assembler
     localIntegralStiffness = std::make_shared<GFE::LocalIntegralStiffness<CompositeBasis,LocalInterpolationRule,RigidBodyMotion> >(bulkCosseratDensity);
   }
 
-  MixedGFEAssembler<CompositeBasis,RigidBodyMotion> mixedAssemblerSmart(compositeBasis,
+  GFE::MixedGFEAssembler<CompositeBasis,RigidBodyMotion> mixedAssemblerSmart(compositeBasis,
                                                                         localIntegralStiffness);
 
   //////////////////////////////////////////////////////////////
diff --git a/test/localprojectedfefunctiontest.cc b/test/localprojectedfefunctiontest.cc
index 2f5f9f5a..9ad47b0b 100644
--- a/test/localprojectedfefunctiontest.cc
+++ b/test/localprojectedfefunctiontest.cc
@@ -66,7 +66,7 @@ evaluateDerivativeFD(const LocalFunction& f, const Dune::FieldVector<ctype, dim>
 
 
 template <int domainDim, int dim>
-void testDerivativeTangentiality(const RealTuple<double,dim>& x,
+void testDerivativeTangentiality(const GFE::RealTuple<double,dim>& x,
                                  const FieldMatrix<double,dim,domainDim>& derivative)
 {
   // By construction, derivatives of RealTuples are always tangent
@@ -74,7 +74,7 @@ void testDerivativeTangentiality(const RealTuple<double,dim>& x,
 
 // the columns of the derivative must be tangential to the manifold
 template <int domainDim, int vectorDim>
-void testDerivativeTangentiality(const UnitVector<double,vectorDim>& x,
+void testDerivativeTangentiality(const GFE::UnitVector<double,vectorDim>& x,
                                  const FieldMatrix<double,vectorDim,domainDim>& derivative)
 {
   for (int i=0; i<domainDim; i++) {
@@ -94,13 +94,13 @@ void testDerivativeTangentiality(const UnitVector<double,vectorDim>& x,
 
 // the columns of the derivative must be tangential to the manifold
 template <int domainDim, int vectorDim>
-void testDerivativeTangentiality(const Rotation<double,vectorDim-1>& x,
+void testDerivativeTangentiality(const GFE::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,typename ... TargetSpaces>
-void testDerivativeTangentiality(const Dune::GFE::ProductManifold<TargetSpaces...>& x,
+void testDerivativeTangentiality(const GFE::ProductManifold<TargetSpaces...>& x,
                                  const FieldMatrix<double,vectorDim,domainDim>& derivative)
 {
   size_t posHelper=0;
@@ -214,7 +214,7 @@ void test(const GeometryType& element)
   std::cout << " --- Testing " << className<TargetSpace>() << ", domain dimension: " << element.dim() << " ---" << std::endl;
 
   std::vector<TargetSpace> testPoints;
-  ValueFactory<TargetSpace>::get(testPoints);
+  GFE::ValueFactory<TargetSpace>::get(testPoints);
 
   int nTestPoints = testPoints.size();
   size_t nVertices = Dune::ReferenceElements<double,domainDim>::general(element).size(domainDim);
@@ -222,7 +222,7 @@ void test(const GeometryType& element)
   // Set up elements of the target space
   std::vector<TargetSpace> corners(nVertices);
 
-  MultiIndex index(nVertices, nTestPoints);
+  GFE::MultiIndex index(nVertices, nTestPoints);
   int numIndices = index.cycle();
 
   for (int i=0; i<numIndices; i++, ++index) {
@@ -260,34 +260,34 @@ int main()
   //  Test functions on 1d elements
   ////////////////////////////////////////////////////////////////
 
-  test<RealTuple<double,1>,1>(GeometryTypes::line);
-  test<UnitVector<double,2>,1>(GeometryTypes::line);
-  test<UnitVector<double,3>,1>(GeometryTypes::line);
-  test<Rotation<double,3>,1>(GeometryTypes::line);
-  typedef Dune::GFE::ProductManifold<RealTuple<double,1>,Rotation<double,3>,UnitVector<double,2> > CrazyManifold;
+  test<GFE::RealTuple<double,1>,1>(GeometryTypes::line);
+  test<GFE::UnitVector<double,2>,1>(GeometryTypes::line);
+  test<GFE::UnitVector<double,3>,1>(GeometryTypes::line);
+  test<GFE::Rotation<double,3>,1>(GeometryTypes::line);
+  typedef GFE::ProductManifold<GFE::RealTuple<double,1>,GFE::Rotation<double,3>,GFE::UnitVector<double,2> > CrazyManifold;
   test<CrazyManifold, 1>(GeometryTypes::line);
 
   ////////////////////////////////////////////////////////////////
   //  Test functions on 2d simplex elements
   ////////////////////////////////////////////////////////////////
 
-  test<RealTuple<double,1>,2>(GeometryTypes::triangle);
-  test<UnitVector<double,2>,2>(GeometryTypes::triangle);
-  test<RealTuple<double,3>,2>(GeometryTypes::triangle);
-  test<UnitVector<double,3>,2>(GeometryTypes::triangle);
-  test<Rotation<double,3>,2>(GeometryTypes::triangle);
-  typedef Dune::GFE::ProductManifold<RealTuple<double,1>,Rotation<double,3>,UnitVector<double,2> > CrazyManifold;
+  test<GFE::RealTuple<double,1>,2>(GeometryTypes::triangle);
+  test<GFE::UnitVector<double,2>,2>(GeometryTypes::triangle);
+  test<GFE::RealTuple<double,3>,2>(GeometryTypes::triangle);
+  test<GFE::UnitVector<double,3>,2>(GeometryTypes::triangle);
+  test<GFE::Rotation<double,3>,2>(GeometryTypes::triangle);
+  typedef GFE::ProductManifold<GFE::RealTuple<double,1>,GFE::Rotation<double,3>,GFE::UnitVector<double,2> > CrazyManifold;
   test<CrazyManifold, 2>(GeometryTypes::triangle);
 
   ////////////////////////////////////////////////////////////////
   //  Test functions on 2d quadrilateral elements
   ////////////////////////////////////////////////////////////////
 
-  test<RealTuple<double,1>,2>(GeometryTypes::quadrilateral);
-  test<UnitVector<double,2>,2>(GeometryTypes::quadrilateral);
-  test<UnitVector<double,3>,2>(GeometryTypes::quadrilateral);
-  test<Rotation<double,3>,2>(GeometryTypes::quadrilateral);
-  typedef Dune::GFE::ProductManifold<RealTuple<double,1>,Rotation<double,3>,UnitVector<double,2> > CrazyManifold;
+  test<GFE::RealTuple<double,1>,2>(GeometryTypes::quadrilateral);
+  test<GFE::UnitVector<double,2>,2>(GeometryTypes::quadrilateral);
+  test<GFE::UnitVector<double,3>,2>(GeometryTypes::quadrilateral);
+  test<GFE::Rotation<double,3>,2>(GeometryTypes::quadrilateral);
+  typedef GFE::ProductManifold<GFE::RealTuple<double,1>,GFE::Rotation<double,3>,GFE::UnitVector<double,2> > CrazyManifold;
   test<CrazyManifold, 2>(GeometryTypes::quadrilateral);
 
 }
diff --git a/test/mixedriemannianpnsolvertest.cc b/test/mixedriemannianpnsolvertest.cc
index 6b297104..cd8c9658 100644
--- a/test/mixedriemannianpnsolvertest.cc
+++ b/test/mixedriemannianpnsolvertest.cc
@@ -56,11 +56,11 @@ using namespace Indices;
 using ValueType = adouble;
 
 //Types for the mixed space
-using DisplacementVector = std::vector<RealTuple<double,dim> >;
-using RotationVector =  std::vector<Rotation<double,dim> >;
+using DisplacementVector = std::vector<GFE::RealTuple<double,dim> >;
+using RotationVector =  std::vector<GFE::Rotation<double,dim> >;
 using Vector = TupleVector<DisplacementVector, RotationVector>;
-using BlockTupleVector = TupleVector<BlockVector<RealTuple<double,dim> >, BlockVector<Rotation<double,dim> > >;
-const int dimRotationTangent = Rotation<double,dim>::TangentVector::dimension;
+using BlockTupleVector = TupleVector<BlockVector<GFE::RealTuple<double,dim> >, BlockVector<GFE::Rotation<double,dim> > >;
+const int dimRotationTangent = GFE::Rotation<double,dim>::TangentVector::dimension;
 
 
 int main (int argc, char *argv[])
@@ -164,18 +164,18 @@ int main (int argc, char *argv[])
                          };
 
   // The target space, with 'double' and 'adouble' as number types
-  using RBM = GFE::ProductManifold<RealTuple<double,dim>,Rotation<double,dim> >;
+  using RBM = GFE::ProductManifold<GFE::RealTuple<double,dim>,GFE::Rotation<double,dim> >;
   using ARBM = typename RBM::template rebind<adouble>::other;
 
   // The total energy
-  auto sumEnergy = std::make_shared<GFE::SumEnergy<CompositeBasis, RealTuple<adouble,dim>,Rotation<adouble,dim> > >();
+  auto sumEnergy = std::make_shared<GFE::SumEnergy<CompositeBasis, GFE::RealTuple<adouble,dim>,GFE::Rotation<adouble,dim> > >();
 
   // The Cosserat shell energy
   using ScalarDeformationLocalFiniteElement = decltype(compositeBasis.localView().tree().child(_0,0).finiteElement());
   using ScalarRotationLocalFiniteElement = decltype(compositeBasis.localView().tree().child(_1,0).finiteElement());
 
-  using AInterpolationRule = std::tuple<LocalGeodesicFEFunction<gridDim, double, ScalarDeformationLocalFiniteElement, RealTuple<adouble,3> >,
-      LocalGeodesicFEFunction<gridDim, double, ScalarRotationLocalFiniteElement, Rotation<adouble,3> > >;
+  using AInterpolationRule = std::tuple<GFE::LocalGeodesicFEFunction<gridDim, double, ScalarDeformationLocalFiniteElement, GFE::RealTuple<adouble,3> >,
+      GFE::LocalGeodesicFEFunction<gridDim, double, ScalarRotationLocalFiniteElement, GFE::Rotation<adouble,3> > >;
 
   auto cosseratDensity = std::make_shared<GFE::PlanarCosseratShellDensity<GridType::Codim<0>::Entity, adouble> >(parameters);
 
@@ -184,14 +184,14 @@ int main (int argc, char *argv[])
   sumEnergy->addLocalEnergy(planarCosseratShellEnergy);
 
   // The Neumann surface load term
-  auto neumannEnergy = std::make_shared<GFE::NeumannEnergy<CompositeBasis, RealTuple<adouble,dim>, Rotation<adouble,dim> > >(neumannBoundary,neumannFunction);
+  auto neumannEnergy = std::make_shared<GFE::NeumannEnergy<CompositeBasis, GFE::RealTuple<adouble,dim>, GFE::Rotation<adouble,dim> > >(neumannBoundary,neumannFunction);
   sumEnergy->addLocalEnergy(neumannEnergy);
 
   // The local assembler
-  LocalGeodesicFEADOLCStiffness<CompositeBasis,RBM> mixedLocalGFEADOLCStiffness(sumEnergy);
+  GFE::LocalGeodesicFEADOLCStiffness<CompositeBasis,RBM> mixedLocalGFEADOLCStiffness(sumEnergy);
 
   // The global assembler
-  MixedGFEAssembler<CompositeBasis, RBM> mixedAssembler(compositeBasis, mixedLocalGFEADOLCStiffness);
+  GFE::MixedGFEAssembler<CompositeBasis, RBM> mixedAssembler(compositeBasis, mixedLocalGFEADOLCStiffness);
 
   using GFEAssemblerWrapper = GFE::GeodesicFEAssemblerWrapper<CompositeBasis, DeformationFEBasis, RBM>;
   GFEAssemblerWrapper assembler(&mixedAssembler, deformationFEBasis);
@@ -239,7 +239,7 @@ int main (int argc, char *argv[])
   const int maxSolverSteps = 1;
   const double initialRegularization = 100;
   const bool instrumented = false;
-  GFE::MixedRiemannianProximalNewtonSolver<CompositeBasis, DeformationFEBasis, RealTuple<double,dim>, DeformationFEBasis, Rotation<double,dim>, BitVector> mixedSolver;
+  GFE::MixedRiemannianProximalNewtonSolver<CompositeBasis, DeformationFEBasis, GFE::RealTuple<double,dim>, DeformationFEBasis, GFE::Rotation<double,dim>, BitVector> mixedSolver;
   mixedSolver.setup(*grid,
                     &mixedAssembler,
                     x,
@@ -252,7 +252,7 @@ int main (int argc, char *argv[])
   mixedSolver.solve();
   x = mixedSolver.getSol();
 
-  RiemannianProximalNewtonSolver<DeformationFEBasis, RBM, GFEAssemblerWrapper> solver;
+  GFE::RiemannianProximalNewtonSolver<DeformationFEBasis, RBM, GFEAssemblerWrapper> solver;
   solver.setup(*grid,
                &assembler,
                xRBM,
diff --git a/test/multiindex.hh b/test/multiindex.hh
index 51b80d4e..960ef39e 100644
--- a/test/multiindex.hh
+++ b/test/multiindex.hh
@@ -3,6 +3,10 @@
 
 #include <vector>
 
+
+namespace Dune::GFE
+{
+
 /** \brief A multi-index
  */
 class MultiIndex
@@ -51,4 +55,6 @@ public:
 
 };
 
+}  // namespace Dune::GFE
+
 #endif
diff --git a/test/nonplanarcosseratenergytest.cc b/test/nonplanarcosseratenergytest.cc
index 6566c170..9be44050 100644
--- a/test/nonplanarcosseratenergytest.cc
+++ b/test/nonplanarcosseratenergytest.cc
@@ -49,7 +49,7 @@ double calculateEnergy(const FlatGridView& flatGridView,
   using namespace Dune::Functions::BasisFactory;
   using namespace Dune::Indices;
 
-  TupleVector<std::vector<RealTuple<double,3> >, std::vector<Rotation<double,3> > > configuration;
+  TupleVector<std::vector<GFE::RealTuple<double,3> >, std::vector<GFE::Rotation<double,3> > > configuration;
   configuration[_0].resize(flatFEBasis.size());
   configuration[_1].resize(flatFEBasis.size());
 
@@ -85,7 +85,7 @@ double calculateEnergy(const FlatGridView& flatGridView,
   auto orientationGridViewFunction = Functions::makeAnalyticGridViewFunction(orientationFunction, curvedGridView);
 
   const auto orientationQuaternionFunction
-    = Functions::makeComposedGridFunction(Rotation<double,3>::matrixToQuaternion,
+    = Functions::makeComposedGridFunction(GFE::Rotation<double,3>::matrixToQuaternion,
                                           orientationGridViewFunction);
 
   BlockVector<FieldVector<double,4> > orientationAsVector(flatFEBasis.size());
@@ -106,7 +106,7 @@ double calculateEnergy(const FlatGridView& flatGridView,
 
   // TODO: Write the curved grid, not the flat one
   // BUG: The second argument should be the displacement, not the deformation
-  CosseratVTKWriter<FlatGridView>::write(flatGridView,
+  GFE::CosseratVTKWriter<FlatGridView>::write(flatGridView,
                                          deformationGridViewFunction,
                                          orientationQuaternionGridViewFunction,
                                          2, // VTK output element order
@@ -119,7 +119,7 @@ double calculateEnergy(const FlatGridView& flatGridView,
   using Element = typename FlatGridView::template Codim<0>::Entity;
   auto density = std::make_shared<GFE::CosseratShellDensity<Element, double> >(materialParameters);
 
-  using ShellEnergy = NonplanarCosseratShellEnergy<FlatFEBasis,
+  using ShellEnergy = GFE::NonplanarCosseratShellEnergy<FlatFEBasis,
       3,                                               // Dimension of the target space
       double,
       GridGeometry>;
@@ -130,7 +130,7 @@ double calculateEnergy(const FlatGridView& flatGridView,
   //  Compute the energy
   ///////////////////////////////////////////////////
 
-  using TargetSpace = GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> >;
+  using TargetSpace = GFE::ProductManifold<GFE::RealTuple<double,3>,GFE::Rotation<double,3> >;
 
   double energy = 0;
 
diff --git a/test/orthogonalmatrixtest.cc b/test/orthogonalmatrixtest.cc
index c1e0ac2c..445cdfd3 100644
--- a/test/orthogonalmatrixtest.cc
+++ b/test/orthogonalmatrixtest.cc
@@ -16,7 +16,7 @@ double eps = 1e-6;
 /** \brief Test whether orthogonal matrix really is orthogonal
  */
 template <class T, int N>
-void testOrthogonality(const OrthogonalMatrix<T,N>& p)
+void testOrthogonality(const GFE::OrthogonalMatrix<T,N>& p)
 {
   Dune::FieldMatrix<T,N,N> prod(0);
   for (int i=0; i<N; i++)
@@ -33,10 +33,10 @@ void testOrthogonality(const OrthogonalMatrix<T,N>& p)
 /** \brief Test orthogonal projection onto the tangent space at p
  */
 template <class T, int N>
-void testProjectionOntoTangentSpace(const OrthogonalMatrix<T,N>& p)
+void testProjectionOntoTangentSpace(const GFE::OrthogonalMatrix<T,N>& p)
 {
   std::vector<FieldMatrix<T,N,N> > testVectors;
-  ValueFactory<FieldMatrix<T,N,N> >::get(testVectors);
+  GFE::ValueFactory<FieldMatrix<T,N,N> >::get(testVectors);
 
   // Test each element in the list
   for (size_t i=0; i<testVectors.size(); i++) {
@@ -70,11 +70,11 @@ void testProjectionOntoTangentSpace(const OrthogonalMatrix<T,N>& p)
 template <class T, int N>
 void test()
 {
-  std::cout << "Testing class " << className<OrthogonalMatrix<T,N> >() << std::endl;
+  std::cout << "Testing class " << className<GFE::OrthogonalMatrix<T,N> >() << std::endl;
 
   // Get set of orthogonal test matrices
-  std::vector<OrthogonalMatrix<T,N> > testPoints;
-  ValueFactory<OrthogonalMatrix<T,N> >::get(testPoints);
+  std::vector<GFE::OrthogonalMatrix<T,N> > testPoints;
+  GFE::ValueFactory<GFE::OrthogonalMatrix<T,N> >::get(testPoints);
 
   int nTestPoints = testPoints.size();
 
diff --git a/test/planarcosseratshelltest.cc b/test/planarcosseratshelltest.cc
index b1a8304c..e323fc4d 100644
--- a/test/planarcosseratshelltest.cc
+++ b/test/planarcosseratshelltest.cc
@@ -43,7 +43,7 @@ int main (int argc, char *argv[])
 {
   MPIHelper::instance(argc, argv);
 
-  using Configuration = TupleVector<std::vector<RealTuple<double,3> >, std::vector<Rotation<double,3> > >;
+  using Configuration = TupleVector<std::vector<GFE::RealTuple<double,3> >, std::vector<GFE::Rotation<double,3> > >;
 
   // solver settings
   const double tolerance                = 1e-4;
@@ -75,7 +75,7 @@ int main (int argc, char *argv[])
   using namespace Dune::Indices;
   using namespace Functions::BasisFactory;
 
-  const int dimRotation = Rotation<double,dim>::embeddedDim;
+  const int dimRotation = GFE::Rotation<double,dim>::embeddedDim;
   auto compositeBasis = makeBasis(
     gridView,
     composite(
@@ -120,8 +120,8 @@ int main (int argc, char *argv[])
     return x;
   });
 
-  BitSetVector<RealTuple<double,dimworld>::TangentVector::dimension> deformationDirichletDofs(deformationFEBasis.size(), false);
-  BitSetVector<Rotation<double,dimworld>::TangentVector::dimension> orientationDirichletDofs(orientationFEBasis.size(), false);
+  BitSetVector<GFE::RealTuple<double,dimworld>::TangentVector::dimension> deformationDirichletDofs(deformationFEBasis.size(), false);
+  BitSetVector<GFE::Rotation<double,dimworld>::TangentVector::dimension> orientationDirichletDofs(orientationFEBasis.size(), false);
 
   const GridView::IndexSet& indexSet = gridView.indexSet();
 
@@ -172,7 +172,7 @@ int main (int argc, char *argv[])
     x[_0][i] = {identity[_0][i][0], identity[_0][i][1], 0.0};
 
   for (auto& rotation : x[_1])
-    rotation = Rotation<double,dimworld>::identity();
+    rotation = GFE::Rotation<double,dimworld>::identity();
 
   //////////////////////////////////
   //  Parameters for the problem
@@ -192,18 +192,18 @@ int main (int argc, char *argv[])
   //////////////////////////////
 
   // The target space, with 'double' and 'adouble' as number types
-  using RigidBodyMotion = GFE::ProductManifold<RealTuple<double,dimworld>,Rotation<double,dimworld> >;
+  using RigidBodyMotion = GFE::ProductManifold<GFE::RealTuple<double,dimworld>,GFE::Rotation<double,dimworld> >;
   using ARigidBodyMotion = typename RigidBodyMotion::template rebind<adouble>::other;
 
   // The total energy
-  auto sumEnergy = std::make_shared<GFE::SumEnergy<CompositeBasis, RealTuple<adouble,dimworld>,Rotation<adouble,dimworld> > >();
+  auto sumEnergy = std::make_shared<GFE::SumEnergy<CompositeBasis, GFE::RealTuple<adouble,dimworld>,GFE::Rotation<adouble,dimworld> > >();
 
   // The Cosserat shell energy
   using ScalarDeformationLocalFiniteElement = decltype(compositeBasis.localView().tree().child(_0,0).finiteElement());
   using ScalarRotationLocalFiniteElement = decltype(compositeBasis.localView().tree().child(_1,0).finiteElement());
 
-  using AInterpolationRule = std::tuple<LocalGeodesicFEFunction<dim, double, ScalarDeformationLocalFiniteElement, RealTuple<adouble,3> >,
-      LocalGeodesicFEFunction<dim, double, ScalarRotationLocalFiniteElement, Rotation<adouble,3> > >;
+  using AInterpolationRule = std::tuple<GFE::LocalGeodesicFEFunction<dim, double, ScalarDeformationLocalFiniteElement, GFE::RealTuple<adouble,3> >,
+      GFE::LocalGeodesicFEFunction<dim, double, ScalarRotationLocalFiniteElement, GFE::Rotation<adouble,3> > >;
 
   auto cosseratDensity = std::make_shared<GFE::PlanarCosseratShellDensity<GridType::Codim<0>::Entity, adouble> >(parameters);
 
@@ -212,17 +212,17 @@ int main (int argc, char *argv[])
   sumEnergy->addLocalEnergy(planarCosseratShellEnergy);
 
   // The Neumann surface load term
-  auto neumannEnergy = std::make_shared<GFE::NeumannEnergy<CompositeBasis, RealTuple<adouble,dimworld>, Rotation<adouble,dimworld> > >(neumannBoundary,neumannFunction);
+  auto neumannEnergy = std::make_shared<GFE::NeumannEnergy<CompositeBasis, GFE::RealTuple<adouble,dimworld>, GFE::Rotation<adouble,dimworld> > >(neumannBoundary,neumannFunction);
   sumEnergy->addLocalEnergy(neumannEnergy);
 
   // The assembler
-  LocalGeodesicFEADOLCStiffness<CompositeBasis,RigidBodyMotion> localGFEADOLCStiffness(sumEnergy);
-  MixedGFEAssembler<CompositeBasis,RigidBodyMotion> mixedAssembler(compositeBasis, localGFEADOLCStiffness);
+  GFE::LocalGeodesicFEADOLCStiffness<CompositeBasis,RigidBodyMotion> localGFEADOLCStiffness(sumEnergy);
+  GFE::MixedGFEAssembler<CompositeBasis,RigidBodyMotion> mixedAssembler(compositeBasis, localGFEADOLCStiffness);
 
-  MixedRiemannianTrustRegionSolver<GridType,
+  GFE::MixedRiemannianTrustRegionSolver<GridType,
       CompositeBasis,
-      DeformationFEBasis, RealTuple<double,dimworld>,
-      OrientationFEBasis, Rotation<double,dimworld> > solver;
+      DeformationFEBasis, GFE::RealTuple<double,dimworld>,
+      OrientationFEBasis, GFE::Rotation<double,dimworld> > solver;
   solver.setup(*grid,
                &mixedAssembler,
                deformationFEBasis,
diff --git a/test/polardecompositiontest.cc b/test/polardecompositiontest.cc
index 7ab8a775..8d8c5618 100644
--- a/test/polardecompositiontest.cc
+++ b/test/polardecompositiontest.cc
@@ -82,7 +82,7 @@ static FieldMatrix<field_type, 3, 3> randomMatrixAlmostOrthogonal(double maxPert
   std::uniform_real_distribution<> dis(0.0, 1.0);   // equally distributed between 0 and upper bound
   for (int i = 0; i < 4; ++i)
     f[i] = dis(gen);
-  Rotation<field_type,3> q(f);
+  GFE::Rotation<field_type,3> q(f);
   q.normalize();
 
   assert(std::abs(1-q.two_norm()) < 1e-12);
diff --git a/test/rotationtest.cc b/test/rotationtest.cc
index f95ac20a..e5b6dedc 100644
--- a/test/rotationtest.cc
+++ b/test/rotationtest.cc
@@ -36,16 +36,16 @@ void testDDExp()
 
         if (j==k) {
 
-          SkewMatrix<double,3> forward(v[i]);
+          GFE::SkewMatrix<double,3> forward(v[i]);
           forward.axial()[j] += eps;
-          Rotation<double,3> forwardQ  = Rotation<double,3>::exp(forward);
+          const auto forwardQ  = GFE::Rotation<double,3>::exp(forward);
 
-          SkewMatrix<double,3> center(v[i]);
-          Rotation<double,3> centerQ   = Rotation<double,3>::exp(center);
+          GFE::SkewMatrix<double,3> center(v[i]);
+          const auto centerQ   = GFE::Rotation<double,3>::exp(center);
 
-          SkewMatrix<double,3> backward(v[i]);
+          GFE::SkewMatrix<double,3> backward(v[i]);
           backward.axial()[j] -= eps;
-          Rotation<double,3> backwardQ = Rotation<double,3>::exp(backward);
+          const auto backwardQ = GFE::Rotation<double,3>::exp(backward);
 
           for (int l=0; l<4; l++)
             fdDDExp[l][j][j] = (forwardQ[l] - 2*centerQ[l] + backwardQ[l]) / (eps*eps);
@@ -53,15 +53,15 @@ void testDDExp()
 
         } else {
 
-          SkewMatrix<double,3> ffV(v[i]);      ffV.axial()[j] += eps;     ffV.axial()[k] += eps;
-          SkewMatrix<double,3> fbV(v[i]);      fbV.axial()[j] += eps;     fbV.axial()[k] -= eps;
-          SkewMatrix<double,3> bfV(v[i]);      bfV.axial()[j] -= eps;     bfV.axial()[k] += eps;
-          SkewMatrix<double,3> bbV(v[i]);      bbV.axial()[j] -= eps;     bbV.axial()[k] -= eps;
+          GFE::SkewMatrix<double,3> ffV(v[i]);      ffV.axial()[j] += eps;     ffV.axial()[k] += eps;
+          GFE::SkewMatrix<double,3> fbV(v[i]);      fbV.axial()[j] += eps;     fbV.axial()[k] -= eps;
+          GFE::SkewMatrix<double,3> bfV(v[i]);      bfV.axial()[j] -= eps;     bfV.axial()[k] += eps;
+          GFE::SkewMatrix<double,3> bbV(v[i]);      bbV.axial()[j] -= eps;     bbV.axial()[k] -= eps;
 
-          Rotation<double,3> forwardForwardQ = Rotation<double,3>::exp(ffV);
-          Rotation<double,3> forwardBackwardQ = Rotation<double,3>::exp(fbV);
-          Rotation<double,3> backwardForwardQ = Rotation<double,3>::exp(bfV);
-          Rotation<double,3> backwardBackwardQ = Rotation<double,3>::exp(bbV);
+          const auto forwardForwardQ = GFE::Rotation<double,3>::exp(ffV);
+          const auto forwardBackwardQ = GFE::Rotation<double,3>::exp(fbV);
+          const auto backwardForwardQ = GFE::Rotation<double,3>::exp(bfV);
+          const auto backwardBackwardQ = GFE::Rotation<double,3>::exp(bbV);
 
           for (int l=0; l<4; l++)
             fdDDExp[l][j][k] = (forwardForwardQ[l] + backwardBackwardQ[l]
@@ -75,7 +75,7 @@ void testDDExp()
 
     // Compute analytical second derivative of exp
     std::array<Dune::FieldMatrix<double,3,3>, 4> ddExp;
-    Rotation<double,3>::DDexp(v[i], ddExp);
+    GFE::Rotation<double,3>::DDexp(v[i], ddExp);
 
     for (int m=0; m<4; m++)
       for (int j=0; j<3; j++)
@@ -89,7 +89,7 @@ void testDDExp()
   }
 }
 
-void testRotation(Rotation<double,3> q)
+void testRotation(GFE::Rotation<double,3> q)
 {
   // Make sure it really is a unit quaternion
   q.normalize();
@@ -111,13 +111,13 @@ void testRotation(Rotation<double,3> q)
   // Turn the matrix back into a quaternion, and check whether it is the same one
   // Since the quaternions form a double covering of SO(3), we may either get q back
   // or -q.  We have to check both.
-  Rotation<double,3> newQ;
+  GFE::Rotation<double,3> newQ;
   newQ.set(matrix);
 
-  Rotation<double,3> diff = newQ;
+  GFE::Rotation<double,3> diff = newQ;
   diff -= q;
 
-  Rotation<double,3> sum  = newQ;
+  GFE::Rotation<double,3> sum  = newQ;
   sum += q;
 
   if (diff.infinity_norm() > 1e-12 && sum.infinity_norm() > 1e-12)
@@ -141,7 +141,7 @@ void testRotation(Rotation<double,3> q)
         for (int l=-2; l<2; l++)
           if (i!=0 || j!=0 || k!=0 || l!=0) {
 
-            Rotation<double,3> q2(Quaternion<double>(i,j,k,l));
+            GFE::Rotation<double,3> q2(GFE::Quaternion<double>(i,j,k,l));
             q2.normalize();
 
             // set up corresponding rotation matrix
@@ -167,7 +167,7 @@ void testRotation(Rotation<double,3> q)
   //   Check the operators 'B' that create an orthonormal basis of H
   // ////////////////////////////////////////////////////////////////
 
-  Quaternion<double> Bq[4];
+  GFE::Quaternion<double> Bq[4];
   Bq[0] = q;
   Bq[1] = q.B(0);
   Bq[2] = q.B(1);
@@ -188,10 +188,10 @@ void testRotation(Rotation<double,3> q)
   //  Check whether the derivativeOfMatrixToQuaternion methods works
   //////////////////////////////////////////////////////////////////////
 
-  Tensor3<double,4,3,3> derivative = Rotation<double,3>::derivativeOfMatrixToQuaternion(matrix);
+  GFE::Tensor3<double,4,3,3> derivative = GFE::Rotation<double,3>::derivativeOfMatrixToQuaternion(matrix);
 
   const double eps = 1e-8;
-  Tensor3<double,4,3,3> derivativeFD;
+  GFE::Tensor3<double,4,3,3> derivativeFD;
 
   for (size_t i=0; i<3; i++)
   {
@@ -202,7 +202,7 @@ void testRotation(Rotation<double,3> q)
       auto backwardMatrix = matrix;
       backwardMatrix[i][j] -= eps;
 
-      Rotation<double,3> forwardRotation, backwardRotation;
+      GFE::Rotation<double,3> forwardRotation, backwardRotation;
       forwardRotation.set(forwardMatrix);
       backwardRotation.set(backwardMatrix);
 
@@ -226,10 +226,10 @@ void testRotation(Rotation<double,3> q)
 }
 
 //! test interpolation between two rotations
-bool testInterpolation(const Rotation<double, 3>& a, const Rotation<double, 3>& b) {
+bool testInterpolation(const GFE::Rotation<double, 3>& a, const GFE::Rotation<double, 3>& b) {
 
   // Compute difference on T_a SO(3)
-  Rotation<double, 3> newB = Rotation<double, 3>::interpolate(a, b, 1.0);
+  auto newB = GFE::Rotation<double, 3>::interpolate(a, b, 1.0);
 
   // Compare matrix representation
   FieldMatrix<double, 3, 3> matB;
@@ -247,8 +247,8 @@ bool testInterpolation(const Rotation<double, 3>& a, const Rotation<double, 3>&
 
 int main (int argc, char *argv[]) try
 {
-  std::vector<Rotation<double,3> > testPoints;
-  ValueFactory<Rotation<double,3> >::get(testPoints);
+  std::vector<GFE::Rotation<double,3> > testPoints;
+  GFE::ValueFactory<GFE::Rotation<double,3> >::get(testPoints);
 
   int nTestPoints = testPoints.size();
 
diff --git a/test/surfacecosseratstressassemblertest.cc b/test/surfacecosseratstressassemblertest.cc
index 429d154f..be51c427 100644
--- a/test/surfacecosseratstressassemblertest.cc
+++ b/test/surfacecosseratstressassemblertest.cc
@@ -61,11 +61,11 @@ static bool sameEntries(FieldVector<double,3> a,FieldVector<double,3> b) {
 }
 
 static bool sameEntries(FieldVector<double,4> a,FieldVector<double,4> b) {
-  Rotation<double,dim> rotationA(a);
-  Rotation<double,dim> rotationB(b);
+  GFE::Rotation<double,dim> rotationA(a);
+  GFE::Rotation<double,dim> rotationB(b);
   rotationA.mult(rotationB);
-  Rotation<double,dim> identity;
-  auto distance = Rotation<double,dim>::distance(rotationA, identity);
+  GFE::Rotation<double,dim> identity;
+  auto distance = GFE::Rotation<double,dim>::distance(rotationA, identity);
   return distance < 0.01;
 }
 
@@ -163,10 +163,10 @@ int main (int argc, char *argv[])
   //               READ IN TEST DATA
   /////////////////////////////////////////////////////////////
 
-  auto deformationMap = Dune::GFE::transformFileToMap<dim>("./stressPlotData/stressPlotTestDeformation");
-  auto initialDeformationMap = Dune::GFE::transformFileToMap<dim>("./stressPlotData/stressPlotTestInitialDeformation");
-  const auto dimRotation = Rotation<double,dim>::embeddedDim;
-  auto rotationMap = Dune::GFE::transformFileToMap<dimRotation>("./stressPlotData/stressPlotTestRotation");
+  auto deformationMap = GFE::transformFileToMap<dim>("./stressPlotData/stressPlotTestDeformation");
+  auto initialDeformationMap = GFE::transformFileToMap<dim>("./stressPlotData/stressPlotTestInitialDeformation");
+  const auto dimRotation = GFE::Rotation<double,dim>::embeddedDim;
+  auto rotationMap = GFE::transformFileToMap<dimRotation>("./stressPlotData/stressPlotTestRotation");
 
   bool deformationIsSymmetric = symmetryTest<dim>(deformationMap, 30);
   bool rotationIsSymmetric = symmetryTest<dimRotation>(rotationMap, 30);
@@ -203,7 +203,7 @@ int main (int argc, char *argv[])
     xInitial[i] += initialDeformationMap.at(stream.str());
   }
 
-  using RotationVector = std::vector<Rotation<double,dim> >;
+  using RotationVector = std::vector<GFE::Rotation<double,dim> >;
   RotationVector rot;
   rot.resize(basisOrderR.size());
   DisplacementVector xOrderR;
@@ -214,7 +214,7 @@ int main (int argc, char *argv[])
   for (std::size_t i = 0; i < basisOrderR.size(); i++) {
     std::stringstream stream;
     stream << xOrderR[i];
-    Rotation<double,dim> rotation(rotationMap.at(stream.str()));
+    GFE::Rotation<double,dim> rotation(rotationMap.at(stream.str()));
     FieldMatrix<double,dim,dim> rotationMatrix(0);
     rotation.matrix(rotationMatrix);
     rot[i].set(rotationMatrix);
@@ -237,7 +237,7 @@ int main (int argc, char *argv[])
   /////////////////////////////////////////////////////////////
 
   auto quadOrder = 4;
-  auto stressAssembler = GFE::SurfaceCosseratStressAssembler<decltype(basisOrderD),decltype(basisOrderR), FieldVector<double,dim>, Rotation<double,dim> >
+  auto stressAssembler = GFE::SurfaceCosseratStressAssembler<decltype(basisOrderD),decltype(basisOrderR), FieldVector<double,dim>, GFE::Rotation<double,dim> >
                            (basisOrderD, basisOrderR);
 
   std::shared_ptr<Elasticity::LocalDensity<dim,ValueType> > elasticDensity;
diff --git a/test/svdtest.cc b/test/svdtest.cc
index e795674a..fb1938d2 100644
--- a/test/svdtest.cc
+++ b/test/svdtest.cc
@@ -16,7 +16,7 @@ int main()
   int nTestValues = 5;
   double maxDiff = 0;
 
-  MultiIndex index(N*M, nTestValues);
+  GFE::MultiIndex index(N*M, nTestValues);
   int numIndices = index.cycle();
 
   for (int i=0; i<numIndices; i++, ++index) {
@@ -32,7 +32,7 @@ int main()
     FieldMatrix<double,N,N> v;
 
     FieldMatrix<double,N,M> testMatrixBackup = testMatrix;
-    svdcmp(testMatrix,w,v);
+    GFE::svdcmp(testMatrix,w,v);
 
     // Multiply the three matrices to see whether we get the original one back
     FieldMatrix<double,N,M> product(0);
diff --git a/test/targetspacetest.cc b/test/targetspacetest.cc
index a0eb6b03..86cafe7c 100644
--- a/test/targetspacetest.cc
+++ b/test/targetspacetest.cc
@@ -231,9 +231,9 @@ void testDerivativeOfHessianOfDistanceSquared(const TargetSpace& a, const Target
   //  Test mixed third derivative with respect to first (once) and second (twice) argument
   /////////////////////////////////////////////////////////////////////////////////////////////
 
-  Tensor3<double,embeddedDim,embeddedDim,embeddedDim> d2d2d2 = TargetSpace::thirdDerivativeOfDistanceSquaredWRTSecondArgument(a, b);
+  GFE::Tensor3<double,embeddedDim,embeddedDim,embeddedDim> d2d2d2 = TargetSpace::thirdDerivativeOfDistanceSquaredWRTSecondArgument(a, b);
 
-  Tensor3<double,embeddedDim,embeddedDim,embeddedDim> d2d2d2_fd;
+  GFE::Tensor3<double,embeddedDim,embeddedDim,embeddedDim> d2d2d2_fd;
 
   for (size_t i=0; i<embeddedDim; i++) {
 
@@ -270,9 +270,9 @@ void testMixedDerivativeOfHessianOfDistanceSquared(const TargetSpace& a, const T
   //  Test mixed third derivative with respect to first (once) and second (twice) argument
   /////////////////////////////////////////////////////////////////////////////////////////////
 
-  Tensor3<double,embeddedDim,embeddedDim,embeddedDim> d1d2d2 = TargetSpace::thirdDerivativeOfDistanceSquaredWRTFirst1AndSecond2Argument(a, b);
+  GFE::Tensor3<double,embeddedDim,embeddedDim,embeddedDim> d1d2d2 = TargetSpace::thirdDerivativeOfDistanceSquaredWRTFirst1AndSecond2Argument(a, b);
 
-  Tensor3<double,embeddedDim,embeddedDim,embeddedDim> d1d2d2_fd;
+  GFE::Tensor3<double,embeddedDim,embeddedDim,embeddedDim> d1d2d2_fd;
 
   for (size_t i=0; i<embeddedDim; i++) {
 
@@ -359,7 +359,7 @@ void test()
   std::cout << "Testing class " << className<TargetSpace>() << std::endl;
 
   std::vector<TargetSpace> testPoints;
-  ValueFactory<TargetSpace>::get(testPoints);
+  GFE::ValueFactory<TargetSpace>::get(testPoints);
 
   int nTestPoints = testPoints.size();
 
@@ -397,23 +397,23 @@ void test()
 int main() try
 {
   // Test the RealTuple class
-  test<RealTuple<double,1> >();
-  test<RealTuple<double,3> >();
+  test<GFE::RealTuple<double,1> >();
+  test<GFE::RealTuple<double,3> >();
 
   // Test the UnitVector class
-  test<UnitVector<double,2> >();
-  test<UnitVector<double,3> >();
-  test<UnitVector<double,4> >();
+  test<GFE::UnitVector<double,2> >();
+  test<GFE::UnitVector<double,3> >();
+  test<GFE::UnitVector<double,4> >();
 
   // Test the rotation class
-  test<Rotation<double,3> >();
+  test<GFE::Rotation<double,3> >();
 
   // Test the ProductManifold class
-  test<Dune::GFE::ProductManifold<RealTuple<double,1>,Rotation<double,3>,UnitVector<double,2> > >();
-  test<Dune::GFE::ProductManifold<Rotation<double,3>,UnitVector<double,5> > >();
+  test<GFE::ProductManifold<GFE::RealTuple<double,1>,GFE::Rotation<double,3>,GFE::UnitVector<double,2> > >();
+  test<GFE::ProductManifold<GFE::Rotation<double,3>,GFE::UnitVector<double,5> > >();
 
   //
-  //     test<HyperbolicHalfspacePoint<double,2> >();
+  //     test<GFE::HyperbolicHalfspacePoint<double,2> >();
 
 }
 catch (Exception& e) {
diff --git a/test/valuefactory.hh b/test/valuefactory.hh
index ed853b9e..2784c0dd 100644
--- a/test/valuefactory.hh
+++ b/test/valuefactory.hh
@@ -10,6 +10,10 @@
 #include <dune/gfe/spaces/rotation.hh>
 #include <dune/gfe/spaces/unitvector.hh>
 
+
+namespace Dune::GFE
+{
+
 /** \brief A class that creates sets of values of various types, to be used in unit tests
  *
  * This is the generic dummy.  The actual work is done in specializations.
@@ -360,5 +364,6 @@ public:
   }
 };
 
+}  // namespace Dune::GFE
 
 #endif
-- 
GitLab