diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 4062ab1f070e2f5a69250f003d193f52b9c67081..3833858e87934e09b964edf05a7e5a0194bcdf62 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -49,6 +49,11 @@ dune_add_test(SOURCES harmonicmaptest.cc
CMAKE_GUARD MPI_FOUND)
endif()
+if (${DUNE_COMMON_VERSION} VERSION_GREATER_EQUAL 2.9)
+ dune_add_test(SOURCES cosseratcontinuumtest.cc
+ TIMEOUT 600)
+endif()
+
dune_add_test(SOURCES geodesicfeassemblerwrappertest.cc
MPI_RANKS 1 4
TIMEOUT 600
diff --git a/test/cosseratcontinuumtest.cc b/test/cosseratcontinuumtest.cc
new file mode 100644
index 0000000000000000000000000000000000000000..0e9aade92703899c67dfc5587a17d02230a1f771
--- /dev/null
+++ b/test/cosseratcontinuumtest.cc
@@ -0,0 +1,255 @@
+#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/parametertree.hh>
+#include <dune/common/typetraits.hh>
+#include <dune/common/bitsetvector.hh>
+
+#include <dune/grid/uggrid.hh>
+#include <dune/grid/utility/structuredgridfactory.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/cosseratenergystiffness.hh>
+#include <dune/gfe/assemblers/mixedlocalgfeadolcstiffness.hh>
+#include <dune/gfe/assemblers/mixedgfeassembler.hh>
+#include <dune/gfe/mixedriemanniantrsolver.hh>
+
+// grid dimension
+const int dim = 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 SolutionType = TupleVector<std::vector<RealTuple<double,dim> >, std::vector<Rotation<double,dim> > >;
+
+ // solver settings
+ const double tolerance = 1e-6;
+ const int maxSolverSteps = 1000;
+ const double initialTrustRegionRadius = 1;
+ const int multigridIterations = 200;
+ const int baseIterations = 100;
+ const double mgTolerance = 1e-10;
+ const double baseTolerance = 1e-8;
+
+ // ///////////////////////////////////////
+ // Create the grid
+ // ///////////////////////////////////////
+ using GridType = UGGrid<dim>;
+
+ //Create a grid with 2 elements, one element will get refined to create different element types, the other one stays a cube
+ std::shared_ptr<GridType> grid = StructuredGridFactory<GridType>::createCubeGrid({0,0,0}, {1,1,1}, {2,1,1});
+
+ // Refine once
+ for (auto&& e : elements(grid->leafGridView())){
+ bool refineThisElement = false;
+ for (int i = 0; i < e.geometry().corners(); i++) {
+ refineThisElement = refineThisElement || (e.geometry().corner(i)[0] > 0.99);
+ }
+ grid->mark(refineThisElement ? 1 : 0, e);
+ }
+ grid->adapt();
+ grid->loadBalance();
+
+ using GridView = GridType::LeafGridView;
+ GridView gridView = grid->leafGridView();
+
+ // ///////////////////////////////////////////
+ // Construct all needed function space bases
+ // ///////////////////////////////////////////
+ using namespace TypeTree::Indices;
+ using namespace Functions::BasisFactory;
+
+ const int dimRotation = Rotation<double,dim>::embeddedDim;
+ auto compositeBasis = makeBasis(
+ gridView,
+ composite(
+ power<dim>(
+ 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 interpolate the identity for both the displacement and the rotation
+ 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<dim> deformationDirichletDofs(deformationFEBasis.size(), false);
+ BitSetVector<dim> 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++) {
+ bool isDirichletDeformation = isDirichlet(identity[_0][i]);
+ for (size_t j=0; j<dim; j++)
+ deformationDirichletDofs[i][j] = isDirichletDeformation;
+ }
+ for (size_t i=0; i<orientationFEBasis.size(); i++) {
+ bool isDirichletOrientation = isDirichlet(identity[_1][i]);
+ for (size_t j=0; j<dim; j++)
+ orientationDirichletDofs[i][j] = isDirichletOrientation;
+ }
+
+ // /////////////////////////////////////////
+ // Determine Neumann dofs and values
+ // /////////////////////////////////////////
+
+ std::function<bool(FieldVector<double,dim>)> isNeumann = [](FieldVector<double,dim> coordinate) {
+ return coordinate[0] > 0.99;
+ };
+
+ 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,dim> values_ = {0,0,0};
+ auto neumannFunction = [&](FieldVector<double, dim>){
+ return values_;
+ };
+
+ // //////////////////////////
+ // Initial iterate
+ // //////////////////////////
+
+ SolutionType x;
+
+ x[_0].resize(compositeBasis.size({0}));
+ x[_1].resize(compositeBasis.size({1}));
+
+ // ////////////////////////////
+ // Parameters for the problem
+ // ////////////////////////////
+
+ ParameterTree parameters;
+ parameters["thickness"] = "0.1";
+ parameters["mu"] = "1";
+ parameters["lambda"] = "1";
+ parameters["mu_c"] = "1";
+ parameters["L_c"] = "0.1";
+ parameters["q"] = "2";
+ parameters["kappa"] = "1";
+ parameters["b1"] = "1";
+ parameters["b2"] = "1";
+ parameters["b3"] = "1";
+
+ // ////////////////////////////
+ // Create an assembler
+ // ////////////////////////////
+
+ CosseratEnergyLocalStiffness<CompositeBasis, dim,adouble> localCosseratEnergy(parameters,
+ neumannBoundary.get(),
+ neumannFunction,
+ nullptr);
+
+ MixedLocalGFEADOLCStiffness<CompositeBasis,
+ RealTuple<double,dim>,
+ Rotation<double,dim> > localGFEADOLCStiffness(&localCosseratEnergy);
+ MixedGFEAssembler<CompositeBasis,
+ RealTuple<double,dim>,
+ Rotation<double,dim> > mixedAssembler(compositeBasis, &localGFEADOLCStiffness);
+
+ MixedRiemannianTrustRegionSolver<GridType,
+ CompositeBasis,
+ DeformationFEBasis, RealTuple<double,dim>,
+ OrientationFEBasis, Rotation<double,dim> > 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.setInitialIterate(x);
+ solver.solve();
+
+ // //////////////////////////////////////////////
+ // Check if solver returns the expected values
+ // //////////////////////////////////////////////
+
+ size_t expectedFinalIteration = 14;
+ double expectedEnergy = 0.67188585;
+
+ 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;
+}