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();
 

test/planarcosseratshelltest.cc

+#include <config.h>
+
+// Includes for the ADOL-C automatic differentiation library
+// Need to come before (almost) all others.
+#include <adolc/adouble.h>
+#include <dune/fufem/utilities/adolcnamespaceinjections.hh>
+
+#include <dune/common/parametertree.hh>
+#include <dune/common/bitsetvector.hh>
+
+#include <dune/grid/io/file/gmshreader.hh>
+#include <dune/grid/uggrid.hh>
+
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <dune/functions/functionspacebases/compositebasis.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
+#include <dune/functions/functionspacebases/interpolate.hh>
+#include <dune/functions/functionspacebases/subspacebasis.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+
+#include <dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh>
+#include <dune/gfe/assemblers/localintegralenergy.hh>
+#include <dune/gfe/assemblers/mixedgfeassembler.hh>
+#include <dune/gfe/assemblers/sumenergy.hh>
+#include <dune/gfe/densities/planarcosseratshelldensity.hh>
+#include <dune/gfe/functions/localgeodesicfefunction.hh>
+#include <dune/gfe/mixedriemanniantrsolver.hh>
+#include <dune/gfe/neumannenergy.hh>
+
+// Dimension of the world space
+const int dimworld = 3;
+
+// Order of the approximation space for the displacement
+const int displacementOrder = 2;
+
+// Order of the approximation space for the microrotations
+const int rotationOrder = 1;
+
+using namespace Dune;
+
+int main (int argc, char *argv[])
+{
+  MPIHelper::instance(argc, argv);
+
+  using Configuration = TupleVector<std::vector<GFE::RealTuple<double,3> >, std::vector<GFE::Rotation<double,3> > >;
+
+  // solver settings
+  const double tolerance                = 1e-4;
+  const int maxSolverSteps              = 20;
+  const double initialTrustRegionRadius = 3.125;
+  const int multigridIterations         = 100;
+  const int baseIterations              = 100;
+  const double mgTolerance              = 1e-10;
+  const double baseTolerance            = 1e-8;
+
+  /////////////////////////////////////////
+  //   Create the grid
+  /////////////////////////////////////////
+
+  const int dim = 2;
+  using GridType = UGGrid<dim>;
+
+  const std::string path = std::string(DUNE_GRID_EXAMPLE_GRIDS_PATH) + "gmsh/";
+  auto grid = GmshReader<GridType>::read(path + "hybrid-testgrid-2d.msh");
+
+  //grid->globalRefine(1);
+
+  using GridView = GridType::LeafGridView;
+  GridView gridView = grid->leafGridView();
+
+  // ///////////////////////////////////////////
+  //  Construct all needed function space bases
+  // ///////////////////////////////////////////
+  using namespace Dune::Indices;
+  using namespace Functions::BasisFactory;
+
+  const int dimRotation = GFE::Rotation<double,dim>::embeddedDim;
+  auto compositeBasis = makeBasis(
+    gridView,
+    composite(
+      power<dimworld>(
+        lagrange<displacementOrder>()
+        ),
+      power<dimRotation>(
+        lagrange<rotationOrder>()
+        )
+      ));
+  using CompositeBasis = decltype(compositeBasis);
+
+  using DeformationFEBasis = Functions::LagrangeBasis<GridView,displacementOrder>;
+  using OrientationFEBasis = Functions::LagrangeBasis<GridView,rotationOrder>;
+
+  DeformationFEBasis deformationFEBasis(gridView);
+  OrientationFEBasis orientationFEBasis(gridView);
+
+  //////////////////////////////////////////////
+  //  Determine Dirichlet dofs
+  //////////////////////////////////////////////
+
+  // This identityBasis is only needed to compute the positions of the Lagrange points
+  auto identityBasis = makeBasis(
+    gridView,
+    composite(
+      power<dim>(
+        lagrange<displacementOrder>()
+        ),
+      power<dim>(
+        lagrange<rotationOrder>()
+        )
+      ));
+  auto deformationPowerBasis = Functions::subspaceBasis(identityBasis, _0);
+  auto rotationPowerBasis = Functions::subspaceBasis(identityBasis, _1);
+
+  MultiTypeBlockVector<BlockVector<FieldVector<double,dim> >,BlockVector<FieldVector<double,dim> > > identity;
+  Functions::interpolate(deformationPowerBasis, identity, [&](FieldVector<double,dim> x){
+    return x;
+  });
+  Functions::interpolate(rotationPowerBasis, identity, [&](FieldVector<double,dim> x){
+    return x;
+  });
+
+  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();
+
+  // Make predicate function that computes which vertices are on the Dirichlet boundary, based on the vertex positions.
+  auto isDirichlet = [](FieldVector<double,dim> coordinate)
+                     {
+                       return coordinate[0] < 0.01;
+                     };
+
+  for (size_t i=0; i<deformationFEBasis.size(); i++)
+    deformationDirichletDofs[i] = isDirichlet(identity[_0][i]);
+
+  for (size_t i=0; i<orientationFEBasis.size(); i++)
+    orientationDirichletDofs[i] = isDirichlet(identity[_1][i]);
+
+  ///////////////////////////////////////////
+  //  Determine Neumann dofs and values
+  ///////////////////////////////////////////
+
+  std::function<bool(FieldVector<double,dim>)> isNeumann
+    = [](FieldVector<double,dim> coordinate)
+      {
+        return coordinate[0] > 0.99 && coordinate[1] > 0.49;
+      };
+
+  BitSetVector<1> neumannVertices(gridView.size(dim), false);
+
+  for (auto&& vertex : vertices(gridView))
+    neumannVertices[indexSet.index(vertex)] = isNeumann(vertex.geometry().corner(0));
+
+  auto neumannBoundary = std::make_shared<BoundaryPatch<GridType::LeafGridView> >(gridView, neumannVertices);
+
+  FieldVector<double,dimworld> values_ = {0,0,60};
+  auto neumannFunction = [&](FieldVector<double, dim>){
+                           return values_;
+                         };
+
+  ////////////////////////////
+  //  Initial iterate
+  ////////////////////////////
+
+  Configuration x;
+
+  x[_0].resize(compositeBasis.size({0}));
+  x[_1].resize(compositeBasis.size({1}));
+
+  for (std::size_t i=0; i<x[_0].size(); ++i)
+    x[_0][i] = {identity[_0][i][0], identity[_0][i][1], 0.0};
+
+  for (auto& rotation : x[_1])
+    rotation = GFE::Rotation<double,dimworld>::identity();
+
+  //////////////////////////////////
+  //  Parameters for the problem
+  //////////////////////////////////
+
+  ParameterTree parameters;
+  parameters["thickness"] = "0.06";
+  parameters["mu"] = "2.7191e+4";
+  parameters["lambda"] = "4.4364e+4";
+  parameters["mu_c"] = "0";
+  parameters["L_c"] = "0.06";
+  parameters["q"] = "2";
+  parameters["kappa"] = "1";         // Shear correction factor
+
+  //////////////////////////////
+  //  Create an assembler
+  //////////////////////////////
+
+  // The target space, with 'double' and 'adouble' as number types
+  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, GFE::RealTuple<adouble,dimworld>,GFE::Rotation<adouble,dimworld> > >();
+
+  // The Cosserat shell energy
+  using AInterpolationRule = std::tuple<GFE::LocalGeodesicFEFunction<DeformationFEBasis, GFE::RealTuple<adouble,3> >,
+      GFE::LocalGeodesicFEFunction<OrientationFEBasis, GFE::Rotation<adouble,3> > >;
+  AInterpolationRule localGFEFunctionA{GFE::LocalGeodesicFEFunction<DeformationFEBasis, GFE::RealTuple<adouble,3> >(deformationFEBasis),
+                                       GFE::LocalGeodesicFEFunction<OrientationFEBasis, GFE::Rotation<adouble,3> >(orientationFEBasis)};
+
+  auto cosseratDensity = std::make_shared<GFE::PlanarCosseratShellDensity<GridType::Codim<0>::Entity, adouble> >(parameters);
+
+  auto planarCosseratShellEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis,AInterpolationRule,ARigidBodyMotion> >(std::move(localGFEFunctionA), cosseratDensity);
+
+  sumEnergy->addLocalEnergy(planarCosseratShellEnergy);
+
+  // The Neumann surface load term
+  auto neumannEnergy = std::make_shared<GFE::NeumannEnergy<CompositeBasis, GFE::RealTuple<adouble,dimworld>, GFE::Rotation<adouble,dimworld> > >(neumannBoundary,neumannFunction);
+  sumEnergy->addLocalEnergy(neumannEnergy);
+
+  // The assembler
+  GFE::LocalGeodesicFEADOLCStiffness<CompositeBasis,RigidBodyMotion> localGFEADOLCStiffness(sumEnergy);
+  GFE::MixedGFEAssembler<CompositeBasis,RigidBodyMotion> mixedAssembler(compositeBasis, localGFEADOLCStiffness);
+
+  GFE::MixedRiemannianTrustRegionSolver<GridType,
+      CompositeBasis,
+      DeformationFEBasis, GFE::RealTuple<double,dimworld>,
+      OrientationFEBasis, GFE::Rotation<double,dimworld> > solver;
+  solver.setup(*grid,
+               &mixedAssembler,
+               deformationFEBasis,
+               orientationFEBasis,
+               x,
+               deformationDirichletDofs,
+               orientationDirichletDofs,
+               tolerance,
+               maxSolverSteps,
+               initialTrustRegionRadius,
+               multigridIterations,
+               mgTolerance,
+               3, 3, 1, // Multigrid V-cycle
+               baseIterations,
+               baseTolerance,
+               false);
+
+  solver.setScaling({1,1,1,0.01,0.01,0.01});
+
+  solver.setInitialIterate(x);
+  solver.solve();
+  x = solver.getSol();
+
+  ////////////////////////////////////////////////
+  //  Check if solver returns the expected values
+  ////////////////////////////////////////////////
+
+  size_t expectedFinalIteration = 9;
+  double expectedEnergy = -6.11748397;
+
+  if (solver.getStatistics().finalIteration != expectedFinalIteration)
+  {
+    std::cerr << "Trust-region solver did " << solver.getStatistics().finalIteration+1
+              << " iterations, instead of the expected '" << expectedFinalIteration+1 << "'!" << std::endl;
+    return 1;
+  }
+
+  if ( std::abs(solver.getStatistics().finalEnergy - expectedEnergy) > 1e-7)
+  {
+    std::cerr << std::setprecision(9);
+    std::cerr << "Final energy is " << solver.getStatistics().finalEnergy
+              << " but '" << expectedEnergy << "' was expected!" << std::endl;
+    return 1;
+  }
+
+  return 0;
+}

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);

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();
 

test/surfacecosseratstressassemblertest.cc

@@ -15,7 +15,11 @@
 #include <dune/common/parametertreeparser.hh>
 #include <dune/common/version.hh>
 
+#if DUNE_VERSION_GTE(DUNE_ELASTICITY, 2, 11)
+#include <dune/elasticity/densities/stvenantkirchhoffdensity.hh>
+#else
 #include <dune/elasticity/materials/stvenantkirchhoffdensity.hh>
+#endif
 
 #include <dune/functions/functionspacebases/interpolate.hh>
 #include <dune/functions/functionspacebases/lagrangebasis.hh>
@@ -27,6 +31,7 @@
 
 #include <dune/gfe/filereader.hh>
 #include <dune/gfe/assemblers/surfacecosseratstressassembler.hh>
+#include <dune/gfe/functions/localgeodesicfefunction.hh>
 #include <dune/gfe/spaces/rotation.hh>
 
 #include <dune/matrix-vector/transpose.hh>
@@ -57,11 +62,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;
 }
 
@@ -159,10 +164,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);
@@ -199,7 +204,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;
@@ -210,7 +215,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);
@@ -233,7 +238,18 @@ int main (int argc, char *argv[])
   /////////////////////////////////////////////////////////////
 
   auto quadOrder = 4;
-  auto stressAssembler = GFE::SurfaceCosseratStressAssembler<decltype(basisOrderD),decltype(basisOrderR), FieldVector<double,dim>, Rotation<double,dim> >
+
+  const Functions::LagrangeBasis<GridView,rotationOrder> scalarBasisR(gridView);
+
+  using LocalGFEFunctionR = GFE::LocalGeodesicFEFunction<decltype(scalarBasisR),GFE::Rotation<double,dim> >;
+
+  LocalGFEFunctionR localGFEFunction(scalarBasisR);
+
+  auto stressAssembler = GFE::SurfaceCosseratStressAssembler<decltype(basisOrderD),
+      decltype(basisOrderR),
+      LocalGFEFunctionR,
+      FieldVector<double,dim>,
+      GFE::Rotation<double,dim> >
                            (basisOrderD, basisOrderR);
 
   std::shared_ptr<Elasticity::LocalDensity<dim,ValueType> > elasticDensity;
@@ -246,7 +262,7 @@ int main (int argc, char *argv[])
   std::vector<FieldMatrix<double,dim,dim> > stressSubstrateCauchyTensor;
   stressAssembler.assembleSubstrateStress<Elasticity::LocalDensity<dim,ValueType> >(x, elasticDensity.get(), quadOrder, stressSubstrate1stPiolaKirchhoffTensor, stressSubstrateCauchyTensor);
   std::vector<FieldMatrix<double,dim,dim> > stressShellBiotTensor;
-  stressAssembler.assembleShellStress(rot, x, xInitial, fLame,/*mu_c*/ 0, surfaceShellBoundary, quadOrder, stressShellBiotTensor);
+  stressAssembler.assembleShellStress(localGFEFunction, rot, x, xInitial, fLame,/*mu_c*/ 0, surfaceShellBoundary, quadOrder, stressShellBiotTensor);
 
   //Now modify ONE value in the rotation function
   int i = 39787;
@@ -256,7 +272,7 @@ int main (int argc, char *argv[])
   Dune::MatrixVector::transpose(transposed, rotationMatrix);
   rot[i].set(transposed);
   std::vector<FieldMatrix<double,dim,dim> > stressShellBiotTensorNotSymmetric;
-  stressAssembler.assembleShellStress(rot, x, xInitial, fLame,/*mu_c*/ 0, surfaceShellBoundary, quadOrder, stressShellBiotTensorNotSymmetric);
+  stressAssembler.assembleShellStress(localGFEFunction, rot, x, xInitial, fLame,/*mu_c*/ 0, surfaceShellBoundary, quadOrder, stressShellBiotTensorNotSymmetric);
   // ... and ONE in the displacement function
   x[i] *= 2;
   std::vector<FieldMatrix<double,dim,dim> > stressSubstrate1stPiolaKirchhoffTensorNotSymmetric;

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);

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) {

test/valuefactory.hh

@@ -10,355 +10,360 @@
 #include <dune/gfe/spaces/rotation.hh>
 #include <dune/gfe/spaces/unitvector.hh>
 
-/** \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.
- */
-template <class T>
-class ValueFactory
+
+namespace Dune::GFE
 {
-public:
-  static void get(std::vector<T>& values);
 
-};
+  /** \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.
+   */
+  template <class T>
+  class ValueFactory
+  {
+  public:
+    static void get(std::vector<T>& values);
 
+  };
 
-/** \brief A class that creates sets of values of various types, to be used in unit tests
- *
- * This is the specialization for RealTuple<1>
- */
-template <>
-class ValueFactory<RealTuple<double,1> >
-{
-public:
-  static void get(std::vector<RealTuple<double,1> >& values) {
 
-    std::vector<double> testPoints = {-3, -1, 0, 2, 4};
+  /** \brief A class that creates sets of values of various types, to be used in unit tests
+   *
+   * This is the specialization for RealTuple<1>
+   */
+  template <>
+  class ValueFactory<RealTuple<double,1> >
+  {
+  public:
+    static void get(std::vector<RealTuple<double,1> >& values) {
 
-    values.resize(testPoints.size());
+      std::vector<double> testPoints = {-3, -1, 0, 2, 4};
 
-    // Set up elements of S^1
-    for (size_t i=0; i<values.size(); i++)
-      values[i] = RealTuple<double,1>(testPoints[i]);
+      values.resize(testPoints.size());
 
-  }
+      // Set up elements of S^1
+      for (size_t i=0; i<values.size(); i++)
+        values[i] = RealTuple<double,1>(testPoints[i]);
 
-};
+    }
 
-/** \brief A class that creates sets of values of various types, to be used in unit tests
- *
- * This is the specialization for RealTuple<3>
- */
-template <>
-class ValueFactory<RealTuple<double,3> >
-{
-public:
-  static void get(std::vector<RealTuple<double,3> >& values) {
+  };
 
-    std::vector<Dune::FieldVector<double,3> > testPoints = {{1,0,0}, {0,1,0}, {-0.838114,0.356751,-0.412667},
-      {-0.490946,-0.306456,0.81551},{-0.944506,0.123687,-0.304319},
-      {-0.6,0.1,-0.2},{0.45,0.12,0.517},
-      {-0.1,0.3,-0.1},{-0.444506,0.123687,0.104319},{-0.7,-0.123687,-0.304319}};
+  /** \brief A class that creates sets of values of various types, to be used in unit tests
+   *
+   * This is the specialization for RealTuple<3>
+   */
+  template <>
+  class ValueFactory<RealTuple<double,3> >
+  {
+  public:
+    static void get(std::vector<RealTuple<double,3> >& values) {
 
-    values.resize(testPoints.size());
+      std::vector<Dune::FieldVector<double,3> > testPoints = {{1,0,0}, {0,1,0}, {-0.838114,0.356751,-0.412667},
+        {-0.490946,-0.306456,0.81551},{-0.944506,0.123687,-0.304319},
+        {-0.6,0.1,-0.2},{0.45,0.12,0.517},
+        {-0.1,0.3,-0.1},{-0.444506,0.123687,0.104319},{-0.7,-0.123687,-0.304319}};
 
-    // Set up elements of S^1
-    for (size_t i=0; i<values.size(); i++)
-      values[i] = RealTuple<double,3>(testPoints[i]);
+      values.resize(testPoints.size());
 
-  }
+      // Set up elements of S^1
+      for (size_t i=0; i<values.size(); i++)
+        values[i] = RealTuple<double,3>(testPoints[i]);
 
-};
+    }
 
-/** \brief A class that creates sets of values of various types, to be used in unit tests
- *
- * This is the specialization for UnitVector<2>
- */
-template <>
-class ValueFactory<UnitVector<double,2> >
-{
-public:
-  static void get(std::vector<UnitVector<double,2> >& values) {
+  };
 
-    std::vector<std::array<double,2> > testPoints = {{1,0}, {0.5,0.5}, {0,1}, {-0.5,0.5}, {-1,0}, {-0.5,-0.5}, {0,-1}, {0.5,-0.5}, {0.1,1}, {1,.1}};
+  /** \brief A class that creates sets of values of various types, to be used in unit tests
+   *
+   * This is the specialization for UnitVector<2>
+   */
+  template <>
+  class ValueFactory<UnitVector<double,2> >
+  {
+  public:
+    static void get(std::vector<UnitVector<double,2> >& values) {
 
-    values.resize(testPoints.size());
+      std::vector<std::array<double,2> > testPoints = {{1,0}, {0.5,0.5}, {0,1}, {-0.5,0.5}, {-1,0}, {-0.5,-0.5}, {0,-1}, {0.5,-0.5}, {0.1,1}, {1,.1}};
 
-    // Set up elements of S^1
-    for (size_t i=0; i<values.size(); i++)
-      values[i] = UnitVector<double,2>(testPoints[i]);
+      values.resize(testPoints.size());
 
-  }
+      // Set up elements of S^1
+      for (size_t i=0; i<values.size(); i++)
+        values[i] = UnitVector<double,2>(testPoints[i]);
 
-};
+    }
 
+  };
 
-/** \brief A class that creates sets of values of various types, to be used in unit tests
- *
- * This is the specialization for UnitVector<3>
- */
-template <>
-class ValueFactory<UnitVector<double,3> >
-{
-public:
-  static void get(std::vector<UnitVector<double,3> >& values) {
 
-    std::vector<std::array<double,3> > testPoints = {{1,0,0}, {0,1,0}, {-0.838114,0.356751,-0.412667},
-      {-0.490946,-0.306456,0.81551},{-0.944506,0.123687,-0.304319},
-      {-0.6,0.1,-0.2},{0.45,0.12,0.517},
-      {-0.1,0.3,-0.1},{-0.444506,0.123687,0.104319},{-0.7,-0.123687,-0.304319}};
+  /** \brief A class that creates sets of values of various types, to be used in unit tests
+   *
+   * This is the specialization for UnitVector<3>
+   */
+  template <>
+  class ValueFactory<UnitVector<double,3> >
+  {
+  public:
+    static void get(std::vector<UnitVector<double,3> >& values) {
 
-    values.resize(testPoints.size());
+      std::vector<std::array<double,3> > testPoints = {{1,0,0}, {0,1,0}, {-0.838114,0.356751,-0.412667},
+        {-0.490946,-0.306456,0.81551},{-0.944506,0.123687,-0.304319},
+        {-0.6,0.1,-0.2},{0.45,0.12,0.517},
+        {-0.1,0.3,-0.1},{-0.444506,0.123687,0.104319},{-0.7,-0.123687,-0.304319}};
 
-    // Set up elements of S^1
-    for (size_t i=0; i<values.size(); i++)
-      values[i] = UnitVector<double,3>(testPoints[i]);
+      values.resize(testPoints.size());
 
-  }
+      // Set up elements of S^1
+      for (size_t i=0; i<values.size(); i++)
+        values[i] = UnitVector<double,3>(testPoints[i]);
 
-};
+    }
 
+  };
 
-/** \brief A class that creates sets of values of various types, to be used in unit tests
- *
- * This is the specialization for UnitVector<4>
- */
-template <>
-class ValueFactory<UnitVector<double,4> >
-{
-public:
-  static void get(std::vector<UnitVector<double,4> >& values) {
 
-    std::vector<std::array<double,4> > testPoints = {{1,0,0,0}, {0,1,0,0}, {-0.838114,0.356751,-0.412667,0.5},
-      {-0.490946,-0.306456,0.81551,0.23},{-0.944506,0.123687,-0.304319,-0.7},
-      {-0.6,0.1,-0.2,0.8},{0.45,0.12,0.517,0},
-      {-0.1,0.3,-0.1,0.73},{-0.444506,0.123687,0.104319,-0.23},{-0.7,-0.123687,-0.304319,0.72}};
+  /** \brief A class that creates sets of values of various types, to be used in unit tests
+   *
+   * This is the specialization for UnitVector<4>
+   */
+  template <>
+  class ValueFactory<UnitVector<double,4> >
+  {
+  public:
+    static void get(std::vector<UnitVector<double,4> >& values) {
 
+      std::vector<std::array<double,4> > testPoints = {{1,0,0,0}, {0,1,0,0}, {-0.838114,0.356751,-0.412667,0.5},
+        {-0.490946,-0.306456,0.81551,0.23},{-0.944506,0.123687,-0.304319,-0.7},
+        {-0.6,0.1,-0.2,0.8},{0.45,0.12,0.517,0},
+        {-0.1,0.3,-0.1,0.73},{-0.444506,0.123687,0.104319,-0.23},{-0.7,-0.123687,-0.304319,0.72}};
 
-    values.resize(testPoints.size());
 
-    // Set up elements of S^1
-    for (size_t i=0; i<values.size(); i++)
-      values[i] = UnitVector<double,4>(testPoints[i]);
+      values.resize(testPoints.size());
 
-  }
+      // Set up elements of S^1
+      for (size_t i=0; i<values.size(); i++)
+        values[i] = UnitVector<double,4>(testPoints[i]);
 
-};
+    }
 
+  };
 
-/** \brief A class that creates sets of values of various types, to be used in unit tests
- *
- * This is the specialization for Rotation<3>
- */
-template <>
-class ValueFactory<Rotation<double,3> >
-{
-public:
-  static void get(std::vector<Rotation<double,3> >& values) {
 
-    std::vector<std::array<double,4> > testPoints = {{1,0,0,0}, {0,1,0,0}, {-0.838114,0.356751,-0.412667,0.5},
-      {-0.490946,-0.306456,0.81551,0.23},{-0.944506,0.123687,-0.304319,-0.7},
-      {-0.6,0.1,-0.2,0.8},{0.45,0.12,0.517,0},
-      {-0.1,0.3,-0.1,0.73},{-0.444506,0.123687,0.104319,-0.23},{-0.7,-0.123687,-0.304319,0.72},
-      {-0.035669, -0.463824, -0.333265, 0.820079}, {0.0178678, 0.916836, 0.367358, 0.155374}};
+  /** \brief A class that creates sets of values of various types, to be used in unit tests
+   *
+   * This is the specialization for Rotation<3>
+   */
+  template <>
+  class ValueFactory<Rotation<double,3> >
+  {
+  public:
+    static void get(std::vector<Rotation<double,3> >& values) {
 
-    values.resize(testPoints.size());
+      std::vector<std::array<double,4> > testPoints = {{1,0,0,0}, {0,1,0,0}, {-0.838114,0.356751,-0.412667,0.5},
+        {-0.490946,-0.306456,0.81551,0.23},{-0.944506,0.123687,-0.304319,-0.7},
+        {-0.6,0.1,-0.2,0.8},{0.45,0.12,0.517,0},
+        {-0.1,0.3,-0.1,0.73},{-0.444506,0.123687,0.104319,-0.23},{-0.7,-0.123687,-0.304319,0.72},
+        {-0.035669, -0.463824, -0.333265, 0.820079}, {0.0178678, 0.916836, 0.367358, 0.155374}};
 
-    // Set up elements of S^1
-    for (std::size_t i=0; i<values.size(); i++)
-      values[i] = Rotation<double,3>(testPoints[i]);
+      values.resize(testPoints.size());
 
-  }
+      // Set up elements of S^1
+      for (std::size_t i=0; i<values.size(); i++)
+        values[i] = Rotation<double,3>(testPoints[i]);
 
-};
+    }
 
-/** \brief A class that creates sets of values of various types, to be used in unit tests
- *
- * This is the specialization for a ProductManifold<RealTuple,Rotation>
- */
-template <>
-class ValueFactory<Dune::GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> > >
-{
-public:
-  static void get(std::vector<Dune::GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> > >& values) {
+  };
 
-    using namespace Dune::Indices;
+  /** \brief A class that creates sets of values of various types, to be used in unit tests
+   *
+   * This is the specialization for a ProductManifold<RealTuple,Rotation>
+   */
+  template <>
+  class ValueFactory<Dune::GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> > >
+  {
+  public:
+    static void get(std::vector<Dune::GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> > >& values) {
 
-    std::vector<RealTuple<double,3> > rValues;
-    ValueFactory<RealTuple<double,3> >::get(rValues);
+      using namespace Dune::Indices;
 
-    std::vector<Rotation<double,3> > qValues;
-    ValueFactory<Rotation<double,3> >::get(qValues);
+      std::vector<RealTuple<double,3> > rValues;
+      ValueFactory<RealTuple<double,3> >::get(rValues);
 
-    int nTestPoints = std::min(rValues.size(), qValues.size());
+      std::vector<Rotation<double,3> > qValues;
+      ValueFactory<Rotation<double,3> >::get(qValues);
 
-    values.resize(nTestPoints);
+      int nTestPoints = std::min(rValues.size(), qValues.size());
 
-    // Set up elements of SE(3)
-    for (int i=0; i<nTestPoints; i++)
-    {
-      values[i][_0] = rValues[i];
-      values[i][_1] = qValues[i];
-    }
+      values.resize(nTestPoints);
 
-  }
+      // Set up elements of SE(3)
+      for (int i=0; i<nTestPoints; i++)
+      {
+        values[i][_0] = rValues[i];
+        values[i][_1] = qValues[i];
+      }
 
-};
+    }
 
-/** \brief A class that creates sets of values of various types, to be used in unit tests
- *
- * This is the specialization for square FieldMatrices
- */
-template <class T, int N>
-class ValueFactory<Dune::FieldMatrix<T,N,N> >
-{
-public:
-  static void get(std::vector<Dune::FieldMatrix<T,N,N> >& values) {
+  };
 
-    int nTestPoints = 10;
-    values.resize(nTestPoints);
+  /** \brief A class that creates sets of values of various types, to be used in unit tests
+   *
+   * This is the specialization for square FieldMatrices
+   */
+  template <class T, int N>
+  class ValueFactory<Dune::FieldMatrix<T,N,N> >
+  {
+  public:
+    static void get(std::vector<Dune::FieldMatrix<T,N,N> >& values) {
 
-    // Set up elements of T^{N \times N}
-    for (int i=0; i<nTestPoints; i++)
-      for (int j=0; j<N; j++)
-        for (int k=0; k<N; k++)
-          values[i][j][k] = std::rand()%100 - 50;
+      int nTestPoints = 10;
+      values.resize(nTestPoints);
 
-  }
+      // Set up elements of T^{N \times N}
+      for (int i=0; i<nTestPoints; i++)
+        for (int j=0; j<N; j++)
+          for (int k=0; k<N; k++)
+            values[i][j][k] = std::rand()%100 - 50;
 
-};
+    }
 
+  };
 
-/** \brief A class that creates sets of values of various types, to be used in unit tests
- *
- * This is the specialization for OrthogonalMatrices
- */
-template <class T, int N>
-class ValueFactory<OrthogonalMatrix<T,N> >
-{
 
-  static Dune::FieldVector<T,N> proj(const Dune::FieldVector<T,N>& u, const Dune::FieldVector<T,N>& v)
+  /** \brief A class that creates sets of values of various types, to be used in unit tests
+   *
+   * This is the specialization for OrthogonalMatrices
+   */
+  template <class T, int N>
+  class ValueFactory<OrthogonalMatrix<T,N> >
   {
-    Dune::FieldVector<T,N> result = u;
-    result *= (v*u) / (u*u);
-    return result;
-  }
 
-public:
-  static void get(std::vector<OrthogonalMatrix<T,N> >& values) {
+    static Dune::FieldVector<T,N> proj(const Dune::FieldVector<T,N>& u, const Dune::FieldVector<T,N>& v)
+    {
+      Dune::FieldVector<T,N> result = u;
+      result *= (v*u) / (u*u);
+      return result;
+    }
+
+  public:
+    static void get(std::vector<OrthogonalMatrix<T,N> >& values) {
+
+      // Get general matrices
+      std::vector<Dune::FieldMatrix<T,N,N> > mValues;
+      ValueFactory<Dune::FieldMatrix<T,N,N> >::get(mValues);
 
-    // Get general matrices
-    std::vector<Dune::FieldMatrix<T,N,N> > mValues;
-    ValueFactory<Dune::FieldMatrix<T,N,N> >::get(mValues);
+      values.resize(mValues.size());
 
-    values.resize(mValues.size());
+      // Do Gram-Schmidt orthogonalization of the rows
+      for (size_t m=0; m<mValues.size(); m++) {
 
-    // Do Gram-Schmidt orthogonalization of the rows
-    for (size_t m=0; m<mValues.size(); m++) {
+        Dune::FieldMatrix<T,N,N>& v = mValues[m];
 
-      Dune::FieldMatrix<T,N,N>& v = mValues[m];
+        if (std::fabs(v.determinant()) < 1e-6)
+          continue;
 
-      if (std::fabs(v.determinant()) < 1e-6)
-        continue;
+        for (int j=0; j<N; j++) {
 
-      for (int j=0; j<N; j++) {
+          for (int i=0; i<j; i++) {
 
-        for (int i=0; i<j; i++) {
+            // v_j = v_j - proj_{v_i} v_j
+            v[j] -= proj(v[i],v[j]);
 
-          // v_j = v_j - proj_{v_i} v_j
-          v[j] -= proj(v[i],v[j]);
+          }
 
+          // normalize
+          v[j] /= v[j].two_norm();
         }
 
-        // normalize
-        v[j] /= v[j].two_norm();
-      }
+        values[m] = OrthogonalMatrix<T,N>(v);
 
-      values[m] = OrthogonalMatrix<T,N>(v);
+      }
 
     }
 
-  }
-
-};
+  };
 
-/** \brief A class that creates sets of values of various types, to be used in unit tests
- *
- * This is the specialization for HyperbolicHalfspacePoint<3>
- */
-template <>
-class ValueFactory<HyperbolicHalfspacePoint<double,2> >
-{
-public:
-  static void get(std::vector<HyperbolicHalfspacePoint<double,2> >& values) {
+  /** \brief A class that creates sets of values of various types, to be used in unit tests
+   *
+   * This is the specialization for HyperbolicHalfspacePoint<3>
+   */
+  template <>
+  class ValueFactory<HyperbolicHalfspacePoint<double,2> >
+  {
+  public:
+    static void get(std::vector<HyperbolicHalfspacePoint<double,2> >& values) {
 
-    std::vector<Dune::FieldVector<double,2> > testPoints = {{0,2}, {0,1}, {0,0.5},
-      {-0.490946,0.81551},{-0.944506,0.304319},
-      {-0.6,0.2},{0.45,0.517},
-      {-0.1,0.1},{-0.444506,0.104319},{-0.7,0.304319}};
+      std::vector<Dune::FieldVector<double,2> > testPoints = {{0,2}, {0,1}, {0,0.5},
+        {-0.490946,0.81551},{-0.944506,0.304319},
+        {-0.6,0.2},{0.45,0.517},
+        {-0.1,0.1},{-0.444506,0.104319},{-0.7,0.304319}};
 
-    values.resize(testPoints.size());
+      values.resize(testPoints.size());
 
-    // Set up elements of S^1
-    for (size_t i=0; i<values.size(); i++)
-      values[i] = HyperbolicHalfspacePoint<double,2>(testPoints[i]);
+      // Set up elements of S^1
+      for (size_t i=0; i<values.size(); i++)
+        values[i] = HyperbolicHalfspacePoint<double,2>(testPoints[i]);
 
-  }
+    }
 
-};
+  };
 
-/** \brief A class that creates sets of values of various types, to be used in unit tests
- *
- * This is the specialization for HyperbolicHalfspacePoint<3>
- */
-template <>
-class ValueFactory<HyperbolicHalfspacePoint<double,3> >
-{
-public:
-  static void get(std::vector<HyperbolicHalfspacePoint<double,3> >& values) {
+  /** \brief A class that creates sets of values of various types, to be used in unit tests
+   *
+   * This is the specialization for HyperbolicHalfspacePoint<3>
+   */
+  template <>
+  class ValueFactory<HyperbolicHalfspacePoint<double,3> >
+  {
+  public:
+    static void get(std::vector<HyperbolicHalfspacePoint<double,3> >& values) {
 
-    std::vector<Dune::FieldVector<double,3> > testPoints = {{1,0,0.01}, {0,1,0.01}, {-0.838114,0.356751,0.412667},
-      {-0.490946,-0.306456,0.81551},{-0.944506,0.123687,0.304319},
-      {-0.6,0.1,0.2},{0.45,0.12,0.517},
-      {-0.1,0.3,0.1},{-0.444506,0.123687,0.104319},{-0.7,-0.123687,0.304319}};
+      std::vector<Dune::FieldVector<double,3> > testPoints = {{1,0,0.01}, {0,1,0.01}, {-0.838114,0.356751,0.412667},
+        {-0.490946,-0.306456,0.81551},{-0.944506,0.123687,0.304319},
+        {-0.6,0.1,0.2},{0.45,0.12,0.517},
+        {-0.1,0.3,0.1},{-0.444506,0.123687,0.104319},{-0.7,-0.123687,0.304319}};
 
-    values.resize(testPoints.size());
+      values.resize(testPoints.size());
 
-    // Set up elements of S^1
-    for (size_t i=0; i<values.size(); i++)
-      values[i] = HyperbolicHalfspacePoint<double,3>(testPoints[i]);
+      // Set up elements of S^1
+      for (size_t i=0; i<values.size(); i++)
+        values[i] = HyperbolicHalfspacePoint<double,3>(testPoints[i]);
 
-  }
+    }
 
-};
+  };
 
 
 
 
 
-/** \brief A class that creates sets of values of various types, to be used in unit tests
- *
- * This is the specialization for ProductManifold<...>
- */
-template <typename ... TargetSpaces>
-class ValueFactory<Dune::GFE::ProductManifold<TargetSpaces...> >
-{
-  using TargetSpace = Dune::GFE::ProductManifold<TargetSpaces...>;
+  /** \brief A class that creates sets of values of various types, to be used in unit tests
+   *
+   * This is the specialization for ProductManifold<...>
+   */
+  template <typename ... TargetSpaces>
+  class ValueFactory<Dune::GFE::ProductManifold<TargetSpaces...> >
+  {
+    using TargetSpace = Dune::GFE::ProductManifold<TargetSpaces...>;
 
-public:
-  static void get(std::vector<TargetSpace >& values) {
+  public:
+    static void get(std::vector<TargetSpace >& values) {
 
-    std::vector<typename TargetSpace::CoordinateType > testPoints(10);
+      std::vector<typename TargetSpace::CoordinateType > testPoints(10);
 
-    std::generate(testPoints.begin(), testPoints.end(), [](){
-      return Dune::GFE::randomFieldVector<typename TargetSpace::field_type,TargetSpace::CoordinateType::dimension>(0.9,1.1) ;
-    });
+      std::generate(testPoints.begin(), testPoints.end(), [](){
+        return Dune::GFE::randomFieldVector<typename TargetSpace::field_type,TargetSpace::CoordinateType::dimension>(0.9,1.1) ;
+      });
 
-    values.resize(testPoints.size());
+      values.resize(testPoints.size());
 
-    std::transform(testPoints.cbegin(),testPoints.cend(),values.begin(),[](const auto& testPoint){return TargetSpace(testPoint);});
-  }
-};
+      std::transform(testPoints.cbegin(),testPoints.cend(),values.begin(),[](const auto& testPoint){return TargetSpace(testPoint);});
+    }
+  };
 
+}  // namespace Dune::GFE
 
 #endif