dune/gfe/sumenergy.hh

+#ifndef DUNE_GFE_SUMENERGY_HH
+#define DUNE_GFE_SUMENERGY_HH
+
+#include <vector>
+
+#include <dune/gfe/localenergy.hh>
+#include <dune/gfe/mixedlocalgeodesicfestiffness.hh>
+
+namespace Dune::GFE {
+
+  /**
+    \brief Assembles the a sum of energies for a single element by summing up the energies of each GFE::LocalEnergy.
+
+           This class works similarly to the class Dune::Elasticity::SumEnergy, where Dune::Elasticity::SumEnergy extends
+           Dune::Elasticity::LocalEnergy and Dune::GFE::SumEnergy extends Dune::GFE::LocalEnergy.
+  */
+
+template<class Basis, class... TargetSpaces>
+class SumEnergy
+  : public Dune::GFE::LocalEnergy<Basis, TargetSpaces...>,
+    public MixedLocalGeodesicFEStiffness<Basis, TargetSpaces...>
+    //Inheriting from MixedLocalGeodesicFEStiffness is hack, and will be replaced eventually; once MixedLocalGFEADOLCStiffness
+    //will be removed and its functionality will be included in LocalGeodesicFEADOLCStiffness this is not needed anymore!
+
+{
+
+  using LocalView = typename Basis::LocalView;
+  using GridView = typename LocalView::GridView;
+  using DT = typename GridView::Grid::ctype;
+  using RT = typename Dune::GFE::LocalEnergy<Basis,TargetSpaces...>::RT;
+
+public:
+  
+  /** \brief Default Constructor*/
+  SumEnergy()
+  {}
+
+  void addLocalEnergy(std::shared_ptr<GFE::LocalEnergy<Basis,TargetSpaces...>> localEnergy)
+  {
+    localEnergies_.push_back(localEnergy);
+  }
+
+  RT energy(const typename Basis::LocalView& localView,
+       const std::vector<TargetSpaces>&... localSolution) const
+  {
+    RT sum = 0.;
+    for ( const auto& localEnergy : localEnergies_ )
+      sum += localEnergy->energy( localView, localSolution... );
+
+    return sum;
+  }
+
+private:
+  // vector containing pointers to the local energies
+  std::vector<std::shared_ptr<GFE::LocalEnergy<Basis, TargetSpaces...>>> localEnergies_;
+};
+
+}  // namespace Dune::GFE
+
+#endif   //#ifndef DUNE_GFE_SUMENERGY_HH

dune/gfe/vtkreader.hh

-// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
-// vi: set et ts=4 sw=2 sts=2:
-#ifndef DUNE_GFE_VTKREADER_HH
-#define DUNE_GFE_VTKREADER_HH
-
-#include <memory>
-
-#include <dune/gfe/vtkfile.hh>
-
-/** \brief Read a grid from a VTK file
- */
-template <class GridType>
-class VTKReader
-{
-public:
-
-  /** \brief Read a grid from a VTK file */
-  static std::unique_ptr<GridType> read(std::string filename)
-  {
-    constexpr auto dimworld = GridType::dimensionworld;
-
-    Dune::GFE::VTKFile vtkFile;
-
-    // Read test file from disk
-    vtkFile.read(filename);
-
-    Dune::GridFactory<GridType> factory;
-
-    for (const auto& v : vtkFile.points_)
-    {
-      // use the first dimworld components as vertex coordinate; discard the rest
-      Dune::FieldVector<typename GridType::ctype, dimworld> pos;
-      for (int i=0; i<dimworld; i++)
-        pos[i] = v[i];
-      factory.insertVertex(pos);
-    }
-
-    for (size_t i=0; i<vtkFile.cellConnectivity_.size(); i+=3)
-    {
-      factory.insertElement(Dune::GeometryTypes::triangle, {vtkFile.cellConnectivity_[i],
-                                       vtkFile.cellConnectivity_[i+1],
-                                       vtkFile.cellConnectivity_[i+2]});
-
-    }
-
-    return std::unique_ptr<GridType>(factory.createGrid());
-  }
-
-};
-
-#endif

src/film-on-substrate.parset → problems/film-on-substrate.parset

@@ -9,13 +9,13 @@ structuredGrid = true
 lower = 0 0 0
 upper = 200 100 200
 
-elements = 10 5 5
+elements = 4 2 2
 
 # Number of grid levels, all elements containing surfaceshell grid vertices will get adaptively refined
-numLevels = 2
+numLevels = 1
 
 # When starting from a file, the stress-free configuration of the surfaceShell is read from a file, this file needs to match the *finest* grid level!
-startFromFile = true
+startFromFile = false
 pathToGridDeformationFile = ./
 gridDeformationFile = deformation
 
@@ -26,7 +26,7 @@ gridDeformation="[1.3*x[0], x[1], x[2]]"
 #  Boundary values
 #############################################
 
-problem = cantilever
+dirichletValues = identity-dirichlet-values
 
 ###  Python predicate specifying all Dirichlet grid vertices
 # x is the vertex coordinate
@@ -34,7 +34,7 @@ dirichletVerticesPredicate = "x[0] < 0.01"
 
 ###  Python predicate specifying all surfaceshell grid vertices, elements conataining these vertices will get adaptively refined
 # x is the vertex coordinate
-surfaceShellVerticesPredicate = "x[2] > 199.99"
+surfaceShellVerticesPredicate = "x[2] > 199.99 and x[0] > 49.99 and x[0] < 150.01"
 
 ###  Python predicate specifying all Neumann grid vertices
 # x is the vertex coordinate
@@ -50,14 +50,31 @@ initialDeformation = "x"
 #  Solver parameters
 #############################################
 
+# Inner solver, cholmod or multigrid
+solvertype = multigrid
+
 # Number of homotopy steps for the Dirichlet boundary conditions
 numHomotopySteps = 1
 
-# Tolerance of the trust region solver
+# Tolerance of the solver
 tolerance = 1e-3
 
-# Max number of steps of the trust region solver
-maxTrustRegionSteps = 20
+# Max number of solver steps
+maxSolverSteps = 1
+
+# Measure convergence
+instrumented = 0
+
+#############################################
+#  Solver parameters specific for proximal newton solver using cholmod
+#############################################
+
+# initial regularization parameter
+initialRegularization = 1000000
+
+#############################################
+#  Solver parameters specific for trust-region solver using multigrid solver 
+#############################################
 
 trustRegionScaling = 1 1 1 0.01 0.01 0.01
 
@@ -85,14 +102,11 @@ mgTolerance = 1e-5
 # Tolerance of the base grid solver
 baseTolerance = 1e-8
 
-# Measure convergence
-instrumented = 0
-
 ############################
 #   Material parameters
 ############################
 
-energy = stvenantkirchhoff
+energy = mooneyrivlin
 
 ## For the Wriggers L-shape example
 [materialParameters]
@@ -108,13 +122,7 @@ surfaceShellParameters = surface-shell-parameters-1-3
 mu_c = 0
 
 # Length scale parameter
-L_c = 0.2 
-
-# Curvature exponent
-q = 2
-
-# Shear correction factor
-kappa = 1
+L_c = 0.2
 
 b1 = 1
 b2 = 1
@@ -140,4 +148,4 @@ mooneyrivlin_energy = log # log, square or ciarlet; different ways to compute th
 # log: Generalized Rivlin model or polynomial hyperelastic model, using  0.5*mooneyrivlin_k*log(det(∇φ)) as the volume-preserving penalty term
 # square: Generalized Rivlin model or polynomial hyperelastic model, using mooneyrivlin_k*(det(∇φ)-1)² as the volume-preserving penalty term
 
-[]
\ No newline at end of file
+[]

src/cantilever-dirichlet-values.py → problems/identity-dirichlet-values.py

@@ -5,13 +5,8 @@ class DirichletValues:
         self.homotopyParameter = homotopyParameter
 
     def deformation(self, x):
-        # Dirichlet b.c. simply clamp the shell in the reference configuration
-        out = x
-
-        return out
-
+        return x
 
     def orientation(self, x):
-        rotation = [[1,0,0], [0, 1, 0], [0, 0, 1]]
+        rotation = [[1,0,0], [0,1,0], [0,0,1]]
         return rotation
-

src/CMakeLists.txt

@@ -5,7 +5,7 @@ set(programs compute-disc-error
              rod3d)
 #            rodobstacle)
 
-if (dune-parmg_FOUND)
+if (dune-parmg_FOUND AND ${DUNE_ELASTICITY_VERSION} VERSION_GREATER_EQUAL 2.8)
 	set(programs film-on-substrate ${programs})
 endif()
 foreach(_program ${programs})

src/cosserat-continuum.cc

@@ -46,7 +46,6 @@
 #include <dune/gfe/nonplanarcosseratshellenergy.hh>
 #include <dune/gfe/cosseratvtkwriter.hh>
 #include <dune/gfe/cosseratvtkreader.hh>
-#include <dune/gfe/vtkreader.hh>
 #include <dune/gfe/geodesicfeassembler.hh>
 #include <dune/gfe/riemanniantrsolver.hh>
 #include <dune/gfe/vertexnormals.hh>
@@ -54,6 +53,10 @@
 #include <dune/gfe/mixedgfeassembler.hh>
 #include <dune/gfe/mixedriemanniantrsolver.hh>
 
+#if HAVE_DUNE_VTK
+#include <dune/vtk/vtkreader.hh>
+#endif
+
 // grid dimension
 const int dim = 2;
 const int dimworld = 2;
@@ -78,45 +81,6 @@ const int blocksize = TargetSpace::TangentVector::dimension;
 
 using namespace Dune;
 
-/** \brief A constant vector-valued function, for simple Neumann boundary values */
-struct NeumannFunction
-    : public Dune::VirtualFunction<FieldVector<double,dimworld>, FieldVector<double,3> >
-{
-    NeumannFunction(const FieldVector<double,3> values,
-                    double homotopyParameter)
-    : values_(values),
-      homotopyParameter_(homotopyParameter)
-    {}
-
-    void evaluate(const FieldVector<double, dimworld>& x, FieldVector<double,3>& out) const {
-        out = 0;
-        out.axpy(homotopyParameter_, values_);
-    }
-
-    FieldVector<double,3> values_;
-    double homotopyParameter_;
-};
-
-/** \brief A constant vector-valued function, for simple volume loads */
-struct VolumeLoad
-    : public Dune::VirtualFunction<FieldVector<double,dimworld>, FieldVector<double,3> >
-{
-    VolumeLoad(const FieldVector<double,3> values,
-               double homotopyParameter)
-    : values_(values),
-      homotopyParameter_(homotopyParameter)
-    {}
-
-    void evaluate(const FieldVector<double, dimworld>& x, FieldVector<double,3>& out) const {
-        out = 0;
-        out.axpy(homotopyParameter_, values_);
-    }
-
-    FieldVector<double,3> values_;
-    double homotopyParameter_;
-};
-
-
 int main (int argc, char *argv[]) try
 {
     // initialize MPI, finalize is done automatically on exit
@@ -203,7 +167,11 @@ int main (int argc, char *argv[]) try
         if (suffix == ".msh")
             grid = std::shared_ptr<GridType>(GmshReader<GridType>::read(path + "/" + gridFile));
         else if (suffix == ".vtu" or suffix == ".vtp")
-            grid = VTKReader<GridType>::read(path + "/" + gridFile);
+#if HAVE_DUNE_VTK
+            grid = VtkReader<GridType>::createGridFromFile(path + "/" + gridFile);
+#else
+            DUNE_THROW(NotImplemented, "Please install dune-vtk for VTK reading support!");
+#endif
     }
 
     grid->globalRefine(numLevels-1);
@@ -218,13 +186,17 @@ int main (int argc, char *argv[]) try
 
     using namespace Dune::Functions::BasisFactory;
 #ifdef MIXED_SPACE
+    const int dimRotation = Rotation<double,dim>::embeddedDim;
     auto compositeBasis = makeBasis(
       gridView,
       composite(
-          lagrange<displacementOrder>(),
-          lagrange<rotationOrder>()
-      )
-    );
+        power<dim>(
+            lagrange<displacementOrder>()
+        ),
+        power<dimRotation>(
+            lagrange<rotationOrder>()
+        )
+    ));
 
     typedef Dune::Functions::LagrangeBasis<GridView,displacementOrder> DeformationFEBasis;
     typedef Dune::Functions::LagrangeBasis<GridView,rotationOrder> OrientationFEBasis;
@@ -427,15 +399,25 @@ int main (int argc, char *argv[]) try
     // ////////////////////////////////////////////////////////////
 
     const ParameterTree& materialParameters = parameterSet.sub("materialParameters");
-    std::shared_ptr<NeumannFunction> neumannFunction;
+    FieldVector<double,3> neumannValues {0,0,0};
     if (parameterSet.hasKey("neumannValues"))
-        neumannFunction = std::make_shared<NeumannFunction>(parameterSet.get<FieldVector<double,3> >("neumannValues"),
-                                                       homotopyParameter);
+        neumannValues = parameterSet.get<FieldVector<double,3> >("neumannValues");
+
+    auto neumannFunction = [&]( FieldVector<double,dim> ) {
+      auto nV = neumannValues;
+      nV *= homotopyParameter;
+      return nV;
+    };
 
-    std::shared_ptr<VolumeLoad> volumeLoad;
+    FieldVector<double,3> volumeLoadValues {0,0,0};
     if (parameterSet.hasKey("volumeLoad"))
-        volumeLoad = std::make_shared<VolumeLoad>(parameterSet.get<FieldVector<double,3> >("volumeLoad"),
-                                                                                          homotopyParameter);
+        volumeLoadValues = parameterSet.get<FieldVector<double,3> >("volumeLoad");
+
+    auto volumeLoad = [&]( FieldVector<double,dim>) {
+      auto vL = volumeLoadValues;
+      vL *= homotopyParameter;
+      return vL;
+    };
 
     if (mpiHelper.rank() == 0) {
         std::cout << "Material parameters:" << std::endl;

test/CMakeLists.txt

@@ -26,5 +26,12 @@ dune_add_test(SOURCES harmonicmaptest.cc
               TIMEOUT 600
               CMAKE_GUARD MPI_FOUND)
 
+if (${DUNE_ELASTICITY_VERSION} VERSION_GREATER_EQUAL 2.7)
+dune_add_test(SOURCES geodesicfeassemblerwrappertest.cc
+              MPI_RANKS 1 4
+              TIMEOUT 600
+              CMAKE_GUARD MPI_FOUND)
+endif()
+
 # Copy the example grid used for testing into the build dir
 file(COPY grids/irregular-square.msh DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/grids)

test/geodesicfeassemblerwrappertest.cc

+#include <config.h>
+
+#include <array>
+
+// 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/typetraits.hh>
+#include <dune/common/bitsetvector.hh>
+
+#include <dune/functions/functionspacebases/compositebasis.hh>
+#include <dune/functions/functionspacebases/interpolate.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
+#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+
+#include <dune/grid/utility/structuredgridfactory.hh>
+#include <dune/grid/uggrid.hh>
+
+#include <dune/gfe/cosseratenergystiffness.hh>
+#include <dune/gfe/geodesicfeassembler.hh>
+#include <dune/gfe/harmonicenergy.hh>
+#include <dune/gfe/localgeodesicfeadolcstiffness.hh>
+#include <dune/gfe/localprojectedfefunction.hh>
+#include <dune/gfe/mixedgfeassembler.hh>
+#include <dune/gfe/mixedlocalgfeadolcstiffness.hh>
+#include <dune/gfe/riemanniantrsolver.hh>
+
+#include <dune/gfe/geodesicfeassemblerwrapper.hh>
+
+#include <dune/istl/multitypeblockmatrix.hh>
+#include <dune/istl/multitypeblockvector.hh>
+
+// grid dimension
+const int gridDim = 2;
+
+// target dimension
+const int dim = 3;
+
+//order of the finite element space
+const int displacementOrder = 2;
+const int rotationOrder = 2;
+
+using namespace Dune;
+using namespace Indices;
+
+//differentiation method: ADOL-C
+using ValueType = adouble;
+
+//Types for the mixed space
+using DisplacementVector = std::vector<RealTuple<double,dim>>;
+using RotationVector =  std::vector<Rotation<double,dim>>;
+using Vector = MultiTypeBlockVector<DisplacementVector, RotationVector>;
+const int dimCR = Rotation<double,dim>::TangentVector::dimension; //dimCorrectionRotation = Dimension of the correction for rotations
+using CorrectionType = MultiTypeBlockVector<BlockVector<FieldVector<double,dim> >, BlockVector<FieldVector<double,dimCR>>>;
+
+using MatrixRow0 = MultiTypeBlockVector<BCRSMatrix<FieldMatrix<double,dim,dim>>,  BCRSMatrix<FieldMatrix<double,dim,dimCR>>>;
+using MatrixRow1 = MultiTypeBlockVector<BCRSMatrix<FieldMatrix<double,dimCR,dim>>, BCRSMatrix<FieldMatrix<double,dimCR,dimCR>>>;
+using MatrixType = MultiTypeBlockMatrix<MatrixRow0,MatrixRow1>;
+
+//Types for the Non-mixed space
+using RBM = RigidBodyMotion<double, dim>;
+const static int blocksize = RBM::TangentVector::dimension;
+using CorrectionTypeWrapped = BlockVector<FieldVector<double, blocksize> >;
+using MatrixTypeWrapped = BCRSMatrix<FieldMatrix<double, blocksize, blocksize> >;
+
+int main (int argc, char *argv[])
+{
+  MPIHelper::instance(argc, argv);
+
+  /////////////////////////////////////////////////////////////////////////
+  //    Create the grid
+  /////////////////////////////////////////////////////////////////////////
+  using GridType = UGGrid<gridDim>;
+  auto grid = StructuredGridFactory<GridType>::createCubeGrid({0,0}, {1,1}, {2,2});
+  grid->globalRefine(2);
+  grid->loadBalance();
+
+  using GridView = GridType::LeafGridView;
+  GridView gridView = grid->leafGridView();
+
+  std::function<bool(FieldVector<double,gridDim>)> isNeumann = [](FieldVector<double,gridDim> coordinate) {
+    return coordinate[0] > 0.99;
+  };
+
+  BitSetVector<1> neumannVertices(gridView.size(gridDim), false);
+  const GridView::IndexSet& indexSet = gridView.indexSet();
+  
+  for (auto&& vertex : vertices(gridView))
+      neumannVertices[indexSet.index(vertex)] = isNeumann(vertex.geometry().corner(0));
+
+  BoundaryPatch<GridView> neumannBoundary(gridView, neumannVertices);
+
+
+  /////////////////////////////////////////////////////////////////////////
+  //  Create a composite basis
+  /////////////////////////////////////////////////////////////////////////
+
+  using namespace Functions::BasisFactory;
+
+  auto compositeBasis = makeBasis(
+    gridView,
+    composite(
+      power<dim>(
+        lagrange<displacementOrder>()
+      ),
+      power<dim>(
+        lagrange<rotationOrder>()
+      )
+  ));
+
+  using CompositeBasis = decltype(compositeBasis);
+  using LocalView = typename CompositeBasis::LocalView;
+
+  /////////////////////////////////////////////////////////////////////////
+  //  Create the energy functions with their parameters
+  /////////////////////////////////////////////////////////////////////////
+
+  //Surface-Cosserat-Energy-Parameters
+  ParameterTree parameters;
+  parameters["thickness"] = "1";
+  parameters["mu"] = "2.7191e+4";
+  parameters["lambda"] = "4.4364e+4";
+  parameters["mu_c"] = "0";
+  parameters["L_c"] = "0.01";
+  parameters["q"] = "2";
+  parameters["kappa"] = "1";
+
+
+  FieldVector<double,dim> values_ = {3e4,2e4,1e4};
+  auto neumannFunction = [&](FieldVector<double, gridDim>){
+    return values_;
+  };
+
+  CosseratEnergyLocalStiffness<decltype(compositeBasis), dim,adouble> cosseratEnergy(parameters,
+                                                                     &neumannBoundary,
+                                                                     neumannFunction,
+                                                                     nullptr);
+  MixedLocalGFEADOLCStiffness<CompositeBasis,
+                              RealTuple<double,dim>,
+                              Rotation<double,dim> > mixedLocalGFEADOLCStiffness(&cosseratEnergy);
+  MixedGFEAssembler<CompositeBasis,
+                    RealTuple<double,dim>,
+                    Rotation<double,dim> > mixedAssembler(compositeBasis, &mixedLocalGFEADOLCStiffness);
+
+  using RBM = RigidBodyMotion<double, dim>;
+  using DeformationFEBasis = Functions::LagrangeBasis<GridView,displacementOrder>;
+  DeformationFEBasis deformationFEBasis(gridView);
+  using GFEAssemblerWrapper = GFE::GeodesicFEAssemblerWrapper<CompositeBasis, DeformationFEBasis, RBM, RealTuple<double, dim>, Rotation<double,dim>>;
+  GFEAssemblerWrapper assembler(&mixedAssembler, deformationFEBasis);
+  
+  /////////////////////////////////////////////////////////////////////////
+  //  Prepare the iterate x where we want to assemble - identity in 2D with z = 0
+  /////////////////////////////////////////////////////////////////////////  
+  auto deformationPowerBasis = makeBasis(
+    gridView,
+    power<gridDim>(
+        lagrange<displacementOrder>()
+  ));
+  BlockVector<FieldVector<double,gridDim> > identity(compositeBasis.size({0}));
+  Functions::interpolate(deformationPowerBasis, identity, [](FieldVector<double,gridDim> x){ return x; });
+  BlockVector<FieldVector<double,dim> > initialDeformation(compositeBasis.size({0}));
+  initialDeformation = 0;
+
+  Vector x;
+  x[_0].resize(compositeBasis.size({0}));
+  x[_1].resize(compositeBasis.size({1}));
+  std::vector<RBM> xRBM(compositeBasis.size({0}));
+  for (int i = 0; i < compositeBasis.size({0}); i++) {
+    for (int j = 0; j < gridDim; j++)
+      initialDeformation[i][j] = identity[i][j];
+    x[_0][i] = initialDeformation[i];
+    for (int j = 0; j < dim; j ++) { // Displacement part
+      xRBM[i].r[j] = x[_0][i][j];
+    }
+    xRBM[i].q = x[_1][i]; // Rotation part
+  }
+
+  //////////////////////////////////////////////////////////////////////////////
+  //  Compute the energy, assemble the Gradient and Hessian using 
+  //  the GeodesicFEAssemblerWrapper and the MixedGFEAssembler and compare!
+  //////////////////////////////////////////////////////////////////////////////
+  CorrectionTypeWrapped gradient;
+  MatrixTypeWrapped hessianMatrix;
+  double energy = assembler.computeEnergy(xRBM);
+  assembler.assembleGradientAndHessian(xRBM, gradient, hessianMatrix, true);
+  double gradientTwoNorm = gradient.two_norm();
+  double gradientInfinityNorm = gradient.infinity_norm();
+  double matrixFrobeniusNorm = hessianMatrix.frobenius_norm();
+
+  CorrectionType gradientMixed;
+  gradientMixed[_0].resize(x[_0].size());
+  gradientMixed[_1].resize(x[_1].size());
+  MatrixType hessianMatrixMixed;
+  double energyMixed = mixedAssembler.computeEnergy(x[_0], x[_1]);
+  mixedAssembler.assembleGradientAndHessian(x[_0], x[_1], gradientMixed[_0], gradientMixed[_1], hessianMatrixMixed, true);
+  double gradientMixedTwoNorm = gradientMixed.two_norm();
+  double gradientMixedInfinityNorm = gradientMixed.infinity_norm();
+  double matrixMixedFrobeniusNorm = hessianMatrixMixed.frobenius_norm();
+
+  if (std::abs(energy - energyMixed)/energyMixed > 1e-8)
+  {
+      std::cerr << std::setprecision(9);
+      std::cerr << "The energy calculated by the GeodesicFEAssemblerWrapper is " << energy << " but "
+                << energyMixed << " (calculated by the MixedGFEAssembler) was expected!" << std::endl;
+      return 1;
+  }
+  if ( std::abs(gradientTwoNorm - gradientMixedTwoNorm)/gradientMixedTwoNorm > 1e-8 ||
+       std::abs(gradientInfinityNorm - gradientMixedInfinityNorm)/gradientMixedInfinityNorm > 1e-8)
+  {
+      std::cerr << std::setprecision(9);
+      std::cerr << "The gradient infinity norm calculated by the GeodesicFEAssemblerWrapper is " << gradientInfinityNorm << " but "
+                << gradientMixedInfinityNorm << " (calculated by the MixedGFEAssembler) was expected!" << std::endl;
+      std::cerr << "The gradient norm calculated by the GeodesicFEAssemblerWrapper is " << gradientTwoNorm << " but "
+                << gradientMixedTwoNorm << " (calculated by the MixedGFEAssembler) was expected!" << std::endl;
+      return 1;
+  }
+
+  if (std::abs(matrixFrobeniusNorm - matrixMixedFrobeniusNorm)/matrixMixedFrobeniusNorm > 1e-8)
+  {
+      std::cerr << std::setprecision(9);
+      std::cerr << "The matrix norm calculated by the GeodesicFEAssemblerWrapper is " << matrixFrobeniusNorm << " but "
+                << matrixMixedFrobeniusNorm << " (calculated by the MixedGFEAssembler) was expected!" << std::endl;
+      return 1;
+  }
+}