diff --git a/src/cosserat-continuum.cc b/src/cosserat-continuum.cc
index 83e5819d975015b486a7b3b6f30d459b02fd3f23..75cb8c3670ced4a94a6e451c8eeba6ce0aed666b 100644
--- a/src/cosserat-continuum.cc
+++ b/src/cosserat-continuum.cc
@@ -501,6 +501,7 @@ int main (int argc, char *argv[]) try
 
             using GFEAssemblerWrapper = Dune::GFE::GeodesicFEAssemblerWrapper<CompositeBasis, DeformationFEBasis, TargetSpace, RealTuple<double, 3>, Rotation<double,3>>;
             GFEAssemblerWrapper assembler(&mixedAssembler, deformationFEBasis);
+            if (parameterSet.get<std::string>("solvertype", "trustRegion") == "trustRegion") {
             RiemannianTrustRegionSolver<DeformationFEBasis, TargetSpace, GFEAssemblerWrapper> solver;
             solver.setup(*grid,
                      &assembler,
@@ -520,6 +521,20 @@ int main (int argc, char *argv[]) try
             solver.setInitialIterate(xTargetSpace);
             solver.solve();
             xTargetSpace = solver.getSol();
+            } else {
+                RiemannianProximalNewtonSolver<DeformationFEBasis, TargetSpace, GFEAssemblerWrapper> solver;
+                solver.setup(*grid,
+                             &assembler,
+                             xTargetSpace,
+                             dirichletDofsTargetSpace,
+                             tolerance,
+                             maxSolverSteps,
+                             initialRegularization,
+                             instrumented);
+                solver.setInitialIterate(xTargetSpace);
+                solver.solve();
+                xTargetSpace = solver.getSol();
+            }
             for (int i = 0; i < xTargetSpace.size(); i++) {
               x[_0][i] = xTargetSpace[i].r;
               x[_1][i] = xTargetSpace[i].q;