diff --git a/dune/gfe/geodesicfeassemblerwrapper.hh b/dune/gfe/geodesicfeassemblerwrapper.hh
new file mode 100644
index 0000000000000000000000000000000000000000..d5b70fd5e58952f5a48659461c5563e40a87f846
--- /dev/null
+++ b/dune/gfe/geodesicfeassemblerwrapper.hh
@@ -0,0 +1,201 @@
+#ifndef GLOBAL_GEODESIC_FE_ASSEMBLER_WRAPPER_HH
+#define GLOBAL_GEODESIC_FE_ASSEMBLER_WRAPPER_HH
+
+#include <dune/gfe/mixedgfeassembler.hh>
+#include <dune/common/tuplevector.hh>
+
+/** \brief A wrapper that wraps a MixedGFEAssembler into an assembler that does not distinguish between the two Finite Element Spaces
+ It reimplements the same methods as these two and transfers the structure of the gradient and the Hessian
+ */
+template <class Basis, class ScalarBasis, class TargetSpace, class MixedSpace0, class MixedSpace1>
+class
+GeodesicFEAssemblerWrapper {
+
+ typedef typename Basis::GridView GridView;
+
+ //! Dimension of the grid.
+ enum { gridDim = GridView::dimension };
+
+ //! Dimension of the tangent space
+ enum { blocksize = TargetSpace::TangentVector::dimension };
+ enum { blocksize0 = MixedSpace0::TangentVector::dimension };
+ enum { blocksize1 = MixedSpace1::TangentVector::dimension };
+
+ //!
+ typedef Dune::FieldMatrix<double, blocksize, blocksize> MatrixBlock;
+ typedef typename MixedGFEAssembler<Basis, MixedSpace0, MixedSpace1>::MatrixType MatrixType;
+
+protected:
+ MixedGFEAssembler<Basis, MixedSpace0, MixedSpace1>* mixedAssembler_;
+
+public:
+ const ScalarBasis& basis_;
+
+ /** \brief Constructor for a given grid */
+ GeodesicFEAssemblerWrapper(MixedGFEAssembler<Basis, MixedSpace0, MixedSpace1>* mixedAssembler, ScalarBasis& basis)
+ : mixedAssembler_(mixedAssembler),
+ basis_(basis)
+ {
+ hessianMixed_ = std::make_unique<MatrixType>();
+ // The two spaces from the mixed version need to have the same embeddedDim as the Target Space
+ assert(MixedSpace0::embeddedDim + MixedSpace1::embeddedDim == TargetSpace::embeddedDim);
+ assert(blocksize0 + blocksize1 == blocksize);
+ }
+
+ /** \brief Assemble the tangent stiffness matrix and the functional gradient together
+ *
+ * This is more efficient than computing them separately, because you need the gradient
+ * anyway to compute the Riemannian Hessian.
+ */
+ virtual void assembleGradientAndHessian(const std::vector<TargetSpace>& sol,
+ Dune::BlockVector<Dune::FieldVector<double, blocksize> >& gradient,
+ Dune::BCRSMatrix<MatrixBlock>& hessian,
+ bool computeOccupationPattern=true) const;
+
+ /** \brief Assemble the gradient */
+ //virtual void assembleGradient(const std::vector<TargetSpace>& sol,
+ // Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const;
+
+ /** \brief Compute the energy of a deformation state */
+ virtual double computeEnergy(const std::vector<TargetSpace>& sol) const;
+
+ /** \brief Get the occupation structure of the Hessian */
+ virtual void getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const;
+
+private:
+ Dune::TupleVector<std::vector<MixedSpace0>,std::vector<MixedSpace1>> splitVector(const std::vector<TargetSpace>& sol) const;
+ std::unique_ptr<MatrixType> hessianMixed_;
+}; // end class
+
+template <class Basis, class ScalarBasis, class TargetSpace, class MixedSpace0, class MixedSpace1>
+Dune::TupleVector<std::vector<MixedSpace0>,std::vector<MixedSpace1>> GeodesicFEAssemblerWrapper<Basis, ScalarBasis, TargetSpace,MixedSpace0,MixedSpace1>::
+splitVector(const std::vector<TargetSpace>& sol) const {
+ using namespace Dune::TypeTree::Indices;
+ // Split the solution into the Deformation and the Rotational part
+ auto n = basis_.size();
+ assert (sol.size() == n);
+ Dune::TupleVector<std::vector<MixedSpace0>,std::vector<MixedSpace1>> solutionSplit;
+ solutionSplit[_0].resize(n);
+ solutionSplit[_1].resize(n);
+ for (int i = 0; i < n; i++) {
+ solutionSplit[_0][i] = sol[i].r; // Deformation part
+ solutionSplit[_1][i] = sol[i].q; // Rotational part
+ }
+ return solutionSplit;
+}
+
+template <class Basis, class ScalarBasis, class TargetSpace, class MixedSpace0, class MixedSpace1>
+void GeodesicFEAssemblerWrapper<Basis, ScalarBasis, TargetSpace,MixedSpace0,MixedSpace1>::
+getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const
+{
+ auto n = basis_.size();
+ nb.resize(n, n);
+ //Retrieve occupation structure for the mixed space and convert it to the non-mixed space
+ //In the mixed space, each index set is for one part of the composite basis
+ //So: nb00 is for the displacement part, nb11 is for the rotation part and nb01 and nb10 (where nb01^T = nb10) is for the coupling
+ Dune::MatrixIndexSet nb00, nb01, nb10, nb11;
+ mixedAssembler_->getMatrixPattern(nb00, nb01, nb10, nb11);
+ auto size0 = nb00.rows();
+ auto size1 = nb11.rows();
+ assert(size0 == size1);
+ assert(size0 == n);
+ //After checking if the sizes match, we can just copy over the occupation pattern
+ //as all occupation patterns work on the same basis combination, so they are all equal
+ nb = nb00;
+}
+
+template <class Basis, class ScalarBasis, class TargetSpace, class MixedSpace0, class MixedSpace1>
+void GeodesicFEAssemblerWrapper<Basis, ScalarBasis, TargetSpace,MixedSpace0,MixedSpace1>::
+assembleGradientAndHessian(const std::vector<TargetSpace>& sol,
+ Dune::BlockVector<Dune::FieldVector<double, blocksize> >& gradient,
+ Dune::BCRSMatrix<MatrixBlock>& hessian,
+ bool computeOccupationPattern) const
+{
+ using namespace Dune::TypeTree::Indices;
+ auto n = basis_.size();
+
+ // Get a split up version of the input
+ auto solutionSplit = splitVector(sol);
+
+ // Define the Matrix and the Gradient in Block Stucture
+ Dune::BlockVector<Dune::FieldVector<double, blocksize0> > gradient0(n);
+ Dune::BlockVector<Dune::FieldVector<double, blocksize1> > gradient1(n);
+
+ if (computeOccupationPattern) {
+ Dune::MatrixIndexSet neighborsPerVertex;
+ getNeighborsPerVertex(neighborsPerVertex);
+ neighborsPerVertex.exportIdx((*hessianMixed_)[_0][_0]);
+ neighborsPerVertex.exportIdx((*hessianMixed_)[_0][_1]);
+ neighborsPerVertex.exportIdx((*hessianMixed_)[_1][_0]);
+ neighborsPerVertex.exportIdx((*hessianMixed_)[_1][_1]);
+ neighborsPerVertex.exportIdx(hessian);
+ }
+
+ *hessianMixed_ = 0;
+ mixedAssembler_->assembleGradientAndHessian(solutionSplit[_0], solutionSplit[_1], gradient0, gradient1, *hessianMixed_, false);
+
+ // Transform gradient and hessian back to the non-mixed structure
+ hessian = 0;
+ gradient.resize(n);
+ gradient = 0;
+ for (int i = 0; i < n; i++) {
+ for (int j = 0; j < blocksize0 + blocksize1; j++)
+ gradient[i][j] = j < blocksize0 ? gradient0[i][j] : gradient1[i][j - blocksize0];
+ }
+
+ // All 4 matrices are nxn;
+ // Each FieldMatrix in (*hessianMixed_)[_0][_0] is blocksize0 x blocksize0
+ // Each FieldMatrix in (*hessianMixed_)[_1][_0] is blocksize1 x blocksize0
+ // Each FieldMatrix in (*hessianMixed_)[_0][_1] is blocksize0 x blocksize1
+ // Each FieldMatrix in (*hessianMixed_)[_1][_1] is blocksize1 x blocksize1
+ // The hessian will then be a nxn matrix where each FieldMatrix is (blocksize0+blocksize1)x(blocksize0+blocksize1)
+
+ for (size_t l = 0; l < blocksize0; l++) {
+ // Extract Upper Left Block of mixed matrix
+ for (auto rowIt = (*hessianMixed_)[_0][_0].begin(); rowIt != (*hessianMixed_)[_0][_0].end(); ++rowIt)
+ for (auto colIt = rowIt->begin(); colIt != rowIt->end(); ++colIt)
+ for(size_t k = 0; k < blocksize0; k++)
+ hessian[rowIt.index()][colIt.index()][k][l] = (*colIt)[k][l];
+ // Extract Lower Left Block of mixed matrix
+ for (auto rowIt = (*hessianMixed_)[_1][_0].begin(); rowIt != (*hessianMixed_)[_1][_0].end(); ++rowIt)
+ for (auto colIt = rowIt->begin(); colIt != rowIt->end(); ++colIt)
+ for(size_t k = 0; k < blocksize1; k++)
+ hessian[rowIt.index()][colIt.index()][k + blocksize0][l] = (*colIt)[k][l];
+ }
+ for (size_t l = 0; l < blocksize1; l++) {
+ // Extract Upper Right Block of mixed matrix
+ for (auto rowIt = (*hessianMixed_)[_0][_1].begin(); rowIt != (*hessianMixed_)[_0][_1].end(); ++rowIt)
+ for (auto colIt = rowIt->begin(); colIt != rowIt->end(); ++colIt)
+ for(size_t k = 0; k < blocksize0; k++)
+ hessian[rowIt.index()][colIt.index()][k][l + blocksize0] = (*colIt)[k][l];
+ // Extract Lower Right Block of mixed matrix
+ for (auto rowIt = (*hessianMixed_)[_1][_1].begin(); rowIt != (*hessianMixed_)[_1][_1].end(); ++rowIt)
+ for (auto colIt = rowIt->begin(); colIt != rowIt->end(); ++colIt)
+ for(size_t k = 0; k < blocksize1; k++)
+ hessian[rowIt.index()][colIt.index()][k + blocksize0][l + blocksize0] = (*colIt)[k][l];
+ }
+}
+
+/*template <class Basis, class ScalarBasis, class TargetSpace, class MixedSpace0, class MixedSpace1>
+void GeodesicFEAssemblerWrapper<Basis, ScalarBasis, TargetSpace,MixedSpace0,MixedSpace1>::
+assembleGradient(const std::vector<TargetSpace>& sol,
+ Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const
+{
+ // Transform gradient back to the non-mixed structure
+ gradient.resize(n);
+ gradient = 0;
+ for (int i = 0; i < n; i++) {
+ for (int j = 0; j < blocksize0 + blocksize1; j++)
+ gradient[i][j] = j < blocksize0 ? gradient0[i][j] : gradient1[i][j - blocksize0];
+ }
+}*/
+
+template <class Basis, class ScalarBasis, class TargetSpace, class MixedSpace0, class MixedSpace1>
+double GeodesicFEAssemblerWrapper<Basis, ScalarBasis, TargetSpace,MixedSpace0,MixedSpace1>::
+computeEnergy(const std::vector<TargetSpace>& sol) const
+{
+ using namespace Dune::TypeTree::Indices;
+ auto solutionSplit = splitVector(sol);
+ return mixedAssembler_->computeEnergy(solutionSplit[_0], solutionSplit[_1]);
+}
+#endif //GLOBAL_GEODESIC_FE_ASSEMBLER_WRAPPER_HH
\ No newline at end of file
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index e28375730ca3df2d1322b4ba0f9d0586ba21583d..7f452954c2ce12749383586c23e009bd6a144542 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -26,5 +26,11 @@ dune_add_test(SOURCES harmonicmaptest.cc
TIMEOUT 600
CMAKE_GUARD MPI_FOUND)
+# Run distributed tests
+dune_add_test(SOURCES geodesicfeassemblerwrappertest.cc
+ MPI_RANKS 1 4
+ TIMEOUT 600
+ CMAKE_GUARD MPI_FOUND)
+
# Copy the example grid used for testing into the build dir
file(COPY grids/irregular-square.msh DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/grids)
diff --git a/test/geodesicfeassemblerwrappertest.cc b/test/geodesicfeassemblerwrappertest.cc
new file mode 100644
index 0000000000000000000000000000000000000000..15a6255ef03fae425d2783eac4f2755043cc39ba
--- /dev/null
+++ b/test/geodesicfeassemblerwrappertest.cc
@@ -0,0 +1,294 @@
+#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/common/parametertree.hh>
+#include <dune/common/parametertreeparser.hh>
+
+#include <dune/elasticity/materials/stvenantkirchhoffdensity.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/fufem/functiontools/boundarydofs.hh>
+
+#include <dune/grid/utility/structuredgridfactory.hh>
+#include <dune/grid/uggrid.hh>
+
+#include <dune/gfe/unitvector.hh>
+#include <dune/gfe/geodesicfeassembler.hh>
+#include <dune/gfe/harmonicenergy.hh>
+#include <dune/gfe/localgeodesicfeadolcstiffness.hh>
+#include <dune/gfe/localgeodesicfefunction.hh>
+#include <dune/gfe/localintegralenergy.hh>
+#include <dune/gfe/localprojectedfefunction.hh>
+#include <dune/gfe/mixedgfeassembler.hh>
+#include <dune/gfe/mixedlocalgfeadolcstiffness.hh>
+#include <dune/gfe/riemanniantrsolver.hh>
+#include <dune/gfe/sumcosseratenergy.hh>
+#include <dune/gfe/surfacecosseratenergy.hh>
+#include <dune/gfe/vertexnormals.hh>
+
+#include <dune/gfe/geodesicfeassemblerwrapper.hh>
+
+#include <dune/istl/multitypeblockmatrix.hh>
+#include <dune/istl/multitypeblockvector.hh>
+
+// grid dimension
+const int dim = 3;
+const int targetDim = 3;
+
+//order of integration
+const int displacementOrder = 2;
+const int rotationOrder = 2;
+
+using namespace Dune;
+using namespace Dune::Indices;
+
+//differentiation method: ADOL-C
+typedef adouble ValueType;
+
+//Types for the mixed space
+typedef std::vector<RealTuple<double,dim> > DisplacementVector;
+typedef std::vector<Rotation<double,dim>> RotationVector;
+typedef MultiTypeBlockVector<DisplacementVector, RotationVector> Vector;
+typedef MultiTypeBlockVector<std::vector<FieldVector<ValueType,dim> >, std::vector<Rotation<ValueType,dim>>> DiffVector; //Equivalent Type used for differentiations
+const int dimCR = Rotation<double,dim>::TangentVector::dimension; //dimCorrectionRotation = Dimension of the correction for rotations
+typedef MultiTypeBlockVector<BlockVector<FieldVector<double,dim> >, BlockVector<FieldVector<double,dimCR>>> CorrectionType;
+
+typedef MultiTypeBlockVector<Dune::BCRSMatrix<Dune::FieldMatrix<double,dim,dim>>, Dune::BCRSMatrix<Dune::FieldMatrix<double,dim,dimCR>>> MatrixRow0;
+typedef MultiTypeBlockVector<Dune::BCRSMatrix<Dune::FieldMatrix<double,dimCR,dim>>, Dune::BCRSMatrix<Dune::FieldMatrix<double,dimCR,dimCR>>> MatrixRow1;
+typedef MultiTypeBlockMatrix<MatrixRow0,MatrixRow1> MatrixType;
+
+//Types for the Non-mixed space
+typedef RigidBodyMotion<double, dim> RBM;
+const static int blocksize = RBM::TangentVector::dimension;
+typedef Dune::BlockVector<Dune::FieldVector<double, blocksize> > CorrectionTypeWrapped;
+typedef Dune::BCRSMatrix<Dune::FieldMatrix<double, blocksize, blocksize> > MatrixTypeWrapped;
+
+
+int main (int argc, char *argv[])
+{
+ MPIHelper::instance(argc, argv);
+
+ FieldVector<double,dim> lower(0), upper(1);
+
+ std::function<bool(Dune::FieldVector<double,dim>)> isSurfaceShell = [](Dune::FieldVector<double,dim> x) {
+ return x[2] > 0.99;
+ };
+
+ /////////////////////////////////////////////////////////////////////////
+ // Create the grid
+ /////////////////////////////////////////////////////////////////////////
+ using GridType = UGGrid<dim>;
+ auto grid = StructuredGridFactory<GridType>::createCubeGrid({0,0,0}, {1,1,1}, {2,1,2});
+ int refine = 1;
+ while(refine > 0) {
+ for (auto&& e : elements(grid->leafGridView())){
+ bool refine = false;
+ for (int i = 0; i < e.geometry().corners(); i++) {
+ refine = refine || isSurfaceShell(e.geometry().corner(i));
+ }
+ grid->mark(refine ? 1 : 0, e);
+ }
+ grid->adapt();
+ refine--;
+ }
+ grid->loadBalance();
+
+ using GridView = GridType::LeafGridView;
+ GridView gridView = grid->leafGridView();
+
+ /////////////////////////////////////////////////////////////////////////
+ // Create a composite basis
+ /////////////////////////////////////////////////////////////////////////
+
+ using namespace Dune::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
+ /////////////////////////////////////////////////////////////////////////
+
+ ParameterTree parameters;
+ //St-Venant-Kirchhoff-Parameters
+ parameters["mu"] = "2.7191e+4";
+ parameters["lambda"] = "4.4364e+4";
+
+ auto elasticDensity = std::make_shared<Elasticity::StVenantKirchhoffDensity<dim,ValueType>>(parameters);
+ auto elasticEnergy = std::make_shared<GFE::LocalIntegralEnergy<LocalView, ValueType>>(elasticDensity);
+
+ //Surface-Cosserat-Energy-Parameters
+ parameters["mu_c"] = "0";
+ parameters["L_c"] = "0.001";
+ parameters["b1"] = "1";
+ parameters["b2"] = "1";
+ parameters["b3"] = "1";
+
+ std::vector<UnitVector<double,dim> > vertexNormals(gridView.size(dim));
+ Dune::FieldVector<double,dim> vertexNormalRaw = {0,0,1};
+ for (int i = 0; i< vertexNormals.size(); i++) {
+ UnitVector vertexNormal(vertexNormalRaw);
+ vertexNormals[i] = vertexNormal;
+ }
+
+ BitSetVector<1> surfaceShellVertices(gridView.size(dim), false);
+ const GridView::IndexSet& indexSet = gridView.indexSet();
+ for (auto&& v : vertices(gridView)) {
+ surfaceShellVertices[indexSet.index(v)] = isSurfaceShell(v.geometry().corner(0));
+ }
+ BoundaryPatch<GridView> surfaceShellBoundary(gridView, surfaceShellVertices);
+
+ typedef MultiLinearGeometry<double, dim-1, dim> ML;
+ std::unordered_map<GridType::GlobalIdSet::IdType, ML> geometriesOnShellBoundary;
+ auto& idSet = grid->globalIdSet();
+
+ for (auto boundaryElement : surfaceShellBoundary) {
+ std::vector<Dune::FieldVector<double,dim>> corners;
+ for (int i = 0; i < boundaryElement.geometry().corners(); i++) {
+ auto corner = boundaryElement.geometry().corner(i);
+ corner[0] *= 1.3;
+ corners.push_back(corner);
+ }
+ GridType::GlobalIdSet::IdType id = idSet.subId(boundaryElement.inside(), boundaryElement.indexInInside(), 1);
+ ML boundaryGeometry(boundaryElement.geometry().type(), corners);
+ geometriesOnShellBoundary.insert({id, boundaryGeometry});
+ }
+
+ std::function<double(Dune::FieldVector<double,dim>)> thicknessF = [](Dune::FieldVector<double,dim> x) {return 1.282;};
+ Dune::FieldVector<double,2> lame(0);
+ lame[0] = 4.22E+06;
+ lame[1] = 8.57E+06;
+ std::function<Dune::FieldVector<double,2>(Dune::FieldVector<double,dim>)> lameF = [lame](Dune::FieldVector<double,dim> x) {return lame;};
+
+ auto surfaceCosseratEnergy = std::make_shared<SurfaceCosseratEnergy<CompositeBasis, DiffVector, ValueType>>(
+ parameters,
+ std::move(vertexNormals),
+ &surfaceShellBoundary,
+ std::move(geometriesOnShellBoundary),
+ thicknessF,
+ lameF);
+ GFE::SumCosseratEnergy<CompositeBasis, DiffVector, ValueType> totalEnergy(elasticEnergy, surfaceCosseratEnergy);
+
+ MixedLocalGFEADOLCStiffness<CompositeBasis,
+ RealTuple<double,dim>,
+ Rotation<double,dim> > mixedLocalGFEADOLCStiffness(&totalEnergy);
+ MixedGFEAssembler<CompositeBasis,
+ RealTuple<double,dim>,
+ Rotation<double,dim> > mixedAssembler(compositeBasis, &mixedLocalGFEADOLCStiffness);
+
+ typedef RigidBodyMotion<double, dim> RBM;
+ typedef Dune::Functions::LagrangeBasis<GridView,displacementOrder> DeformationFEBasis;
+ DeformationFEBasis deformationFEBasis(gridView);
+ typedef GeodesicFEAssemblerWrapper<CompositeBasis, DeformationFEBasis, RBM, RealTuple<double, dim>, Rotation<double,dim>> GFEAssemblerWrapper;
+ GFEAssemblerWrapper assembler(&mixedAssembler, deformationFEBasis);
+
+ /////////////////////////////////////////////////////////////////////////
+ // Prepare the iterate x where we want to assemble
+ /////////////////////////////////////////////////////////////////////////
+ auto deformationPowerBasis = makeBasis(
+ gridView,
+ power<dim>(
+ lagrange<displacementOrder>()
+ ));
+ BlockVector<FieldVector<double,dim> > identity(compositeBasis.size({0}));
+ Dune::Functions::interpolate(deformationPowerBasis, identity, [](FieldVector<double,dim> x){ return x; });
+
+ 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++) {
+ x[_0][i] = identity[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 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();
+
+ //These values were taken from assembling using a fufem-basis
+ double expectedEnergy = 563243.074;
+ double expectedGradientTwoNorm = 2857597.54;
+ double expectedGradientInfinityNorm = 1251651.27;
+ double expectedMatrixFrobeniusNorm = 1.63129339e+09;
+
+ if ( std::abs(energy - expectedEnergy)/expectedEnergy > 1e-4 || std::abs(energy - energyMixed)/energyMixed > 1e-4)
+ {
+ std::cerr << std::setprecision(9);
+ std::cerr << "The energy calculated by the GeodesicFEAssemblerWrapper is " << energy << " but "
+ << expectedEnergy << " was expected!" << std::endl;
+ std::cerr << "The energy calculated by the MixedGFEAssembler is " << energyMixed << "!" << std::endl;
+ return 1;
+ }
+ if ( std::abs(gradientTwoNorm - expectedGradientTwoNorm)/expectedGradientTwoNorm > 1e-4 ||
+ std::abs(gradientInfinityNorm - expectedGradientInfinityNorm)/expectedGradientInfinityNorm > 1e-8 ||
+ std::abs(gradientTwoNorm - gradientMixedTwoNorm)/gradientMixedTwoNorm > 1e-4 ||
+ 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 "
+ << expectedGradientInfinityNorm << " was expected!" << std::endl;
+ std::cerr << "The gradient infinity norm calculated by the MixedGFEAssembler is " << gradientMixedInfinityNorm << "!" << std::endl;
+ std::cerr << "The gradient norm calculated by the GeodesicFEAssemblerWrapper is " << gradientTwoNorm << " but "
+ << expectedGradientTwoNorm << " was expected!" << std::endl;
+ std::cerr << "The gradient norm calculated by the MixedGFEAssembler is " << gradientMixedTwoNorm << "!" << std::endl;
+ return 1;
+ }
+
+ if ( std::abs(matrixFrobeniusNorm - expectedMatrixFrobeniusNorm)/expectedMatrixFrobeniusNorm > 1e-4 ||
+ std::abs(matrixFrobeniusNorm - matrixMixedFrobeniusNorm)/matrixMixedFrobeniusNorm > 1e-4 )
+ {
+ std::cerr << std::setprecision(9);
+ std::cerr << "The matrix norm calculated by the GeodesicFEAssemblerWrapper is " << matrixFrobeniusNorm << " but "
+ << expectedMatrixFrobeniusNorm << " was expected!" << std::endl;
+ std::cerr << "The matrix norm calculated by the MixedGFEAssembler is " << matrixMixedFrobeniusNorm << "!" << std::endl;
+ return 1;
+ }
+}