diff --git a/dune/gfe/mixedriemannianpnsolver.cc b/dune/gfe/mixedriemannianpnsolver.cc
new file mode 100644
index 0000000000000000000000000000000000000000..eea8d8cbe295aa3186f6d59bdb223e00f7dcbd75
--- /dev/null
+++ b/dune/gfe/mixedriemannianpnsolver.cc
@@ -0,0 +1,381 @@
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/timer.hh>
+
+#include <dune/gfe/parallel/matrixcommunicator.hh>
+#include <dune/gfe/parallel/vectorcommunicator.hh>
+
+template <class MixedBasis,
+ class Basis0,
+ class TargetSpace0,
+ class Basis1,
+ class TargetSpace1,
+ class BitVector>
+void Dune::GFE::MixedRiemannianProximalNewtonSolver<MixedBasis,Basis0,TargetSpace0,Basis1,TargetSpace1,BitVector>::
+setup(const GridType& grid,
+ const MixedGFEAssembler<MixedBasis, TargetSpace0, TargetSpace1>* assembler,
+ const SolutionType& x,
+ const BitVector& dirichletNodes,
+ double tolerance,
+ int maxProximalNewtonSteps,
+ double initialRegularization,
+ bool instrumented)
+{
+ grid_ = &grid;
+ assembler_ = assembler;
+ x_ = x;
+ tolerance_ = tolerance;
+ maxProximalNewtonSteps_ = maxProximalNewtonSteps;
+ initialRegularization_ = initialRegularization;
+
+ //////////////////////////////////////////////////////////////////
+ // Create global numbering for matrix and vector transfer
+ //////////////////////////////////////////////////////////////////
+#if HAVE_MPI
+ globalMapper0_ = std::make_unique<GlobalMapper0>(grid_->leafGridView());
+ globalMapper1_ = std::make_unique<GlobalMapper1>(grid_->leafGridView());
+
+ // Transfer all Dirichlet data to the master processor
+ using VectorCommunicator0 = VectorCommunicator<GlobalMapper0, typename GridType::LeafGridView::Communication, Dune::BitSetVector<blocksize0> >;
+ VectorCommunicator0 vectorComm0(*globalMapper0_,
+ grid_->leafGridView().comm(),
+ 0);
+ using VectorCommunicator1 = VectorCommunicator<GlobalMapper1, typename GridType::LeafGridView::Communication, Dune::BitSetVector<blocksize1> >;
+ VectorCommunicator1 vectorComm1(*globalMapper1_,
+ grid_->leafGridView().comm(),
+ 0);
+
+ using namespace Dune::Indices;
+ Dune::BitSetVector<blocksize0> dirichletNodes0(dirichletNodes[_0].size(),false);
+ Dune::BitSetVector<blocksize1> dirichletNodes1(dirichletNodes[_1].size(),false);
+ for (std::size_t i = 0; i < dirichletNodes[_0].size(); i++)
+ for (int j = 0; j < blocksize0; j++)
+ dirichletNodes0[i][j] = dirichletNodes[_0][i][j];
+ for (std::size_t i = 0; i < dirichletNodes[_1].size(); i++)
+ for (int j = 0; j < blocksize1; j++)
+ dirichletNodes1[i][j] = dirichletNodes[_1][i][j];
+
+ auto globalDirichletNodes0 = vectorComm0.reduceCopy(dirichletNodes0);
+ auto globalDirichletNodes1 = vectorComm1.reduceCopy(dirichletNodes1);
+
+ BitVector globalDirichletNodes01;
+ globalDirichletNodes01[_0].resize(globalDirichletNodes0.size());
+ globalDirichletNodes01[_1].resize(globalDirichletNodes1.size());
+ for (std::size_t i = 0; i < globalDirichletNodes0.size(); i++)
+ for (int j = 0; j < blocksize0; j++)
+ globalDirichletNodes01[_0][i][j] = globalDirichletNodes0[i][j];
+ for (std::size_t i = 0; i < globalDirichletNodes1.size(); i++)
+ for (int j = 0; j < blocksize1; j++)
+ globalDirichletNodes01[_1][i][j] = globalDirichletNodes1[i][j];
+
+ auto globalDirichletNodes = new BitVector(globalDirichletNodes01);
+#else
+ auto globalDirichletNodes = new BitVector(dirichletNodes);
+#endif
+
+ //////////////////////////////////////////////////////////////////
+ // Create the inner solver using a direct solver
+ //////////////////////////////////////////////////////////////////
+
+ innerSolver_ = std::make_shared<Solvers::CholmodSolver<MatrixType,CorrectionType,BitVector> >();
+
+ innerSolver_->setIgnore(*globalDirichletNodes);
+ hessianMatrix_ = std::make_unique<MatrixType>();
+}
+
+
+template <class MixedBasis,
+ class Basis0,
+ class TargetSpace0,
+ class Basis1,
+ class TargetSpace1,
+ class BitVector>
+void Dune::GFE::MixedRiemannianProximalNewtonSolver<MixedBasis,Basis0,TargetSpace0,Basis1,TargetSpace1,BitVector>::solve()
+{
+ int rank = grid_->comm().rank();
+
+ // /////////////////////////////////////////////////////
+ // Proximal Newton Solver
+ // /////////////////////////////////////////////////////
+
+ using namespace Dune::TypeTree::Indices;
+
+ Dune::Timer energyTimer;
+ double oldEnergy = assembler_->computeEnergy(x_[_0], x_[_1]);
+ if (this->verbosity_ == Solver::FULL)
+ std::cout << "Energy computation took " << energyTimer.elapsed() << " sec." << std::endl;
+
+ oldEnergy = grid_->comm().sum(oldEnergy);
+
+ bool recomputeGradientHessian = true;
+ CorrectionType rhs, rhs_global;
+ MatrixType stiffnessMatrix;
+#if HAVE_MPI
+ using VectorCommunicator0 = VectorCommunicator<GlobalMapper0, typename GridType::LeafGridView::Communication, CorrectionType0>;
+ VectorCommunicator0 vectorComm0(*globalMapper0_,
+ grid_->leafGridView().comm(),
+ 0);
+ using VectorCommunicator1 = VectorCommunicator<GlobalMapper1, typename GridType::LeafGridView::Communication, CorrectionType1>;
+ VectorCommunicator1 vectorComm1(*globalMapper1_,
+ grid_->leafGridView().comm(),
+ 0);
+
+ LocalMapper0 localMapper0 = MapperFactory<Basis0>::createLocalMapper(grid_->leafGridView());
+ LocalMapper1 localMapper1 = MapperFactory<Basis1>::createLocalMapper(grid_->leafGridView());
+
+ MatrixCommunicator<GlobalMapper0,
+ typename GridType::LeafGridView,
+ typename GridType::LeafGridView,
+ MatrixType00,
+ LocalMapper0,
+ LocalMapper0> matrixComm00(*globalMapper0_,
+ grid_->leafGridView(),
+ localMapper0,
+ localMapper0,
+ 0);
+ MatrixCommunicator<GlobalMapper1,
+ typename GridType::LeafGridView,
+ typename GridType::LeafGridView,
+ MatrixType11,
+ LocalMapper1,
+ LocalMapper1> matrixComm11(*globalMapper1_,
+ grid_->leafGridView(),
+ localMapper1,
+ localMapper1,
+ 0);
+ MatrixCommunicator<GlobalMapper0,
+ typename GridType::LeafGridView,
+ typename GridType::LeafGridView,
+ MatrixType01,
+ LocalMapper0,
+ LocalMapper1,
+ GlobalMapper1> matrixComm01(*globalMapper0_,
+ *globalMapper1_,
+ grid_->leafGridView(),
+ grid_->leafGridView(),
+ localMapper0,
+ localMapper1,
+ 0);
+ MatrixCommunicator<GlobalMapper1,
+ typename GridType::LeafGridView,
+ typename GridType::LeafGridView,
+ MatrixType10,
+ LocalMapper1,
+ LocalMapper0,
+ GlobalMapper0> matrixComm10(*globalMapper1_,
+ *globalMapper0_,
+ grid_->leafGridView(),
+ grid_->leafGridView(),
+ localMapper1,
+ localMapper0,
+ 0);
+#endif
+ double totalAssemblyTime = 0.0;
+ double totalSolverTime = 0.0;
+ double regularization = initialRegularization_;
+ for (std::size_t i=0; i<maxProximalNewtonSteps_; i++)
+ {
+
+ Dune::Timer totalTimer;
+ if (this->verbosity_ == Solver::FULL and rank==0) {
+ std::cout << "----------------------------------------------------" << std::endl;
+ std::cout << " Mixed Proximal Newton Step Number: " << i
+ << ", regularization parameter: " << regularization
+ << ", energy: " << oldEnergy << std::endl;
+ std::cout << "----------------------------------------------------" << std::endl;
+ }
+
+ CorrectionType corr;
+ corr[_0].resize(x_[_0].size());
+ corr[_1].resize(x_[_1].size());
+ corr = 0;
+
+ if (recomputeGradientHessian) {
+ Dune::Timer assemblyTimer;
+ assembler_->assembleGradientAndHessian(x_[_0],
+ x_[_1],
+ rhs[_0],
+ rhs[_1],
+ *hessianMatrix_,
+ i==0 // assemble occupation pattern only for the first call
+ );
+
+ rhs *= -1; // The right hand side is the _negative_ gradient
+
+ if (this->verbosity_ == Solver::FULL)
+ std::cout << "Assembly took " << assemblyTimer.elapsed() << " sec." << std::endl;
+ totalAssemblyTime += assemblyTimer.elapsed();
+
+#if HAVE_MPI
+ // Transfer matrix data
+ stiffnessMatrix[_0][_0] = matrixComm00.reduceAdd((*hessianMatrix_)[_0][_0]);
+ stiffnessMatrix[_0][_1] = matrixComm01.reduceAdd((*hessianMatrix_)[_0][_1]);
+ stiffnessMatrix[_1][_0] = matrixComm10.reduceAdd((*hessianMatrix_)[_1][_0]);
+ stiffnessMatrix[_1][_1] = matrixComm11.reduceAdd((*hessianMatrix_)[_1][_1]);
+
+ // Transfer vector data
+ rhs_global[_0] = vectorComm0.reduceAdd(rhs[_0]);
+ rhs_global[_1] = vectorComm1.reduceAdd(rhs[_1]);
+#else
+ stiffnessMatrix = *hessianMatrix_;
+ rhs_global = rhs;
+#endif
+ recomputeGradientHessian = false;
+ }
+
+ CorrectionType corr_global;
+ corr_global[_0].resize(rhs_global[_0].size());
+ corr_global[_1].resize(rhs_global[_1].size());
+ corr_global = 0;
+ bool solvedByInnerSolver = true;
+
+ if (rank==0)
+ {
+ // Add the regularization - Identity Matrix for now
+ for (std::size_t i=0; i<stiffnessMatrix[_0][_0].N(); i++)
+ for (int j=0; j<blocksize0; j++)
+ stiffnessMatrix[_0][_0][i][i][j][j] += regularization;
+ for (std::size_t i=0; i<stiffnessMatrix[_1][_1].N(); i++)
+ for (int j=0; j<blocksize1; j++)
+ stiffnessMatrix[_1][_1][i][i][j][j] += regularization;
+
+ innerSolver_->setProblem(stiffnessMatrix,corr_global,rhs_global);
+ innerSolver_->preprocess();
+
+ ///////////////////////////////
+ // Solve !
+ ///////////////////////////////
+
+ std::cout << "Solve quadratic problem using cholmod solver..." << std::endl;
+ Dune::Timer solutionTimer;
+ try {
+ innerSolver_->solve();
+ } catch (Dune::Exception &e) {
+ std::cerr << "Error while solving: " << e << std::endl;
+ solvedByInnerSolver = false;
+ corr_global = 0;
+ }
+ std::cout << "Solving the quadratic problem took " << solutionTimer.elapsed() << " seconds." << std::endl;
+ totalSolverTime += solutionTimer.elapsed();
+ }
+#if HAVE_MPI
+ // Distribute solution
+ if (grid_->comm().size()>1 and rank==0)
+ std::cout << "Transfer solution back to root process ..." << std::endl;
+
+ corr[_0] = vectorComm0.scatter(corr_global[_0]);
+ corr[_1] = vectorComm1.scatter(corr_global[_1]);
+#else
+ corr = corr_global;
+#endif
+ double corrNorm = corr.infinity_norm();
+ double corrGlobalInfinityNorm = grid_->comm().max(corrNorm);
+ if (std::isnan(corrGlobalInfinityNorm))
+ solvedByInnerSolver = false;
+
+ double energy = 0;
+ double modelDecrease = 0;
+ SolutionType newIterate = x_;
+ if (i == maxProximalNewtonSteps_ - 1)
+ std::cout << i+1 << " proximal newton steps were taken, the maximum was reached." << std::endl << "Total solver time: " << totalSolverTime << " sec., total assembly time: " << totalAssemblyTime << " sec." << std::endl;
+
+ if (solvedByInnerSolver) {
+ if (this->verbosity_ == NumProc::FULL && rank==0)
+ std::cout << "Infinity norm of the correction: " << corrGlobalInfinityNorm << std::endl;
+
+ if (corrGlobalInfinityNorm < tolerance_) {
+ if (verbosity_ == NumProc::FULL and rank==0)
+ std::cout << "CORRECTION IS SMALL ENOUGH" << std::endl;
+
+ if (verbosity_ != NumProc::QUIET and rank==0)
+ std::cout << i+1 << " proximal newton steps were taken" << std::endl << "Total solver time: " << totalSolverTime << " sec., total assembly time: " << totalAssemblyTime << " sec." << std::endl;
+ break;
+ }
+
+ // ////////////////////////////////////////////////////
+ // Check whether proximal newton step can be accepted
+ // ////////////////////////////////////////////////////
+
+ for (size_t j=0; j<newIterate[_0].size(); j++)
+ newIterate[_0][j] = TargetSpace0::exp(newIterate[_0][j], corr[_0][j]);
+
+ for (size_t j=0; j<newIterate[_1].size(); j++)
+ newIterate[_1][j] = TargetSpace1::exp(newIterate[_1][j], corr[_1][j]);
+ try {
+ energy = assembler_->computeEnergy(newIterate[_0],newIterate[_1]);
+ } catch (Dune::Exception &e) {
+ std::cerr << "Error while computing the energy of the new Iterate: " << e << std::endl;
+ std::cerr << "Redoing proximal newton step with higher regularization parameter ..." << std::endl;
+ solvedByInnerSolver = false;
+ }
+ solvedByInnerSolver = grid_->comm().min(solvedByInnerSolver);
+ if (!solvedByInnerSolver) {
+ newIterate = x_;
+ energy = oldEnergy;
+ } else {
+ energy = grid_->comm().sum(energy);
+
+ // compute the model decrease
+ // It is $ m(x) - m(x+s) = -<g,s> - 0.5 <s, Hs>
+ // Note that rhs = -g
+ CorrectionType tmp(corr);
+ hessianMatrix_->mv(corr,tmp);
+ modelDecrease = rhs*corr - 0.5 * (corr*tmp);
+ modelDecrease = grid_->comm().sum(modelDecrease);
+
+ double relativeModelDecrease = modelDecrease / std::fabs(energy);
+
+ if (verbosity_ == NumProc::FULL and rank==0) {
+ std::cout << "Absolute model decrease: " << modelDecrease
+ << ", functional decrease: " << oldEnergy - energy << std::endl;
+ std::cout << "Relative model decrease: " << relativeModelDecrease
+ << ", functional decrease: " << (oldEnergy - energy)/energy << std::endl;
+ }
+
+ assert(modelDecrease >= 0);
+
+ if (energy >= oldEnergy and rank==0) {
+ if (this->verbosity_ == NumProc::FULL)
+ std::cout << "Direction is not a descent direction!" << std::endl;
+ }
+ }
+ }
+ // //////////////////////////////////////////////
+ // Check for acceptance of the step
+ // //////////////////////////////////////////////
+ if ( solvedByInnerSolver && oldEnergy >= energy && (oldEnergy-energy) / modelDecrease > 0.9)
+ {
+ // very successful iteration
+
+ x_ = newIterate;
+ regularization *= 0.5;
+
+ // current energy becomes 'oldEnergy' for the next iteration
+ oldEnergy = energy;
+
+ recomputeGradientHessian = true;
+
+ } else if ((solvedByInnerSolver && oldEnergy >= energy && (oldEnergy-energy) / modelDecrease > 0.01)
+ || std::abs(oldEnergy-energy) < 1e-12) {
+ // successful iteration
+ x_ = newIterate;
+
+ // current energy becomes 'oldEnergy' for the next iteration
+ oldEnergy = energy;
+
+ recomputeGradientHessian = true;
+
+ } else {
+
+ // unsuccessful iteration
+
+ // Increase the regularization parameter
+ regularization *= 2;
+
+ if (this->verbosity_ == NumProc::FULL and rank==0)
+ std::cout << "Unsuccessful iteration!" << std::endl;
+ }
+
+ if (rank==0)
+ std::cout << "iteration took " << totalTimer.elapsed() << " sec." << std::endl;
+ }
+}
diff --git a/dune/gfe/mixedriemannianpnsolver.hh b/dune/gfe/mixedriemannianpnsolver.hh
new file mode 100644
index 0000000000000000000000000000000000000000..a0314bbd02f5dc22db8807d9e7a3f299f5cfe608
--- /dev/null
+++ b/dune/gfe/mixedriemannianpnsolver.hh
@@ -0,0 +1,118 @@
+#ifndef DUNE_GFE_MIXED_RIEMANNIAN_PROXIMAL_NEWTON_SOLVER_HH
+#define DUNE_GFE_MIXED_RIEMANNIAN_PROXIMAL_NEWTON_SOLVER_HH
+
+#include <vector>
+
+#include <dune/common/bitsetvector.hh>
+
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/bvector.hh>
+#include <dune/istl/multitypeblockmatrix.hh>
+
+#include <dune/grid/utility/globalindexset.hh>
+
+#include <dune/solvers/solvers/cholmodsolver.hh>
+
+#include <dune/gfe/parallel/mapperfactory.hh>
+
+
+namespace Dune::GFE
+{
+
+ /** \brief Riemannian proximal Newton solver for geodesic finite-element problems */
+ template <class MixedBasis,
+ class Basis0,
+ class TargetSpace0,
+ class Basis1,
+ class TargetSpace1,
+ class BitVector>
+ class MixedRiemannianProximalNewtonSolver
+ : public NumProc
+ {
+ using GridType = typename MixedBasis::GridView::Grid;
+
+ const static int blocksize0 = TargetSpace0::TangentVector::dimension;
+ const static int blocksize1 = TargetSpace1::TangentVector::dimension;
+
+ const static int gridDim = GridType::dimension;
+
+ // Centralize the field type here
+ using field_type = double;
+
+ using MatrixType00 = BCRSMatrix<FieldMatrix<field_type, blocksize0, blocksize0> >;
+ using MatrixType01 = BCRSMatrix<FieldMatrix<field_type, blocksize0, blocksize1> >;
+ using MatrixType10 = BCRSMatrix<FieldMatrix<field_type, blocksize1, blocksize0> >;
+ using MatrixType11 = BCRSMatrix<FieldMatrix<field_type, blocksize1, blocksize1> >;
+ typedef MultiTypeBlockMatrix<MultiTypeBlockVector<MatrixType00,MatrixType01>,
+ MultiTypeBlockVector<MatrixType10,MatrixType11> > MatrixType;
+
+ using CorrectionType0 = BlockVector<FieldVector<field_type, blocksize0> >;
+ using CorrectionType1 = BlockVector<FieldVector<field_type, blocksize1> >;
+ using CorrectionType = MultiTypeBlockVector<CorrectionType0, CorrectionType1>;
+ using SolutionType = TupleVector<std::vector<TargetSpace0>, std::vector<TargetSpace1> >;
+#if HAVE_MPI
+ typedef typename MapperFactory<Basis0>::GlobalMapper GlobalMapper0;
+ typedef typename MapperFactory<Basis1>::GlobalMapper GlobalMapper1;
+ typedef typename MapperFactory<Basis0>::LocalMapper LocalMapper0;
+ typedef typename MapperFactory<Basis1>::LocalMapper LocalMapper1;
+#endif
+
+ public:
+
+ MixedRiemannianProximalNewtonSolver()
+ : NumProc(NumProc::FULL)
+ {}
+
+ void setup(const GridType& grid,
+ const MixedGFEAssembler<MixedBasis, TargetSpace0, TargetSpace1>* assembler,
+ const SolutionType& x,
+ const BitVector& dirichletNodes,
+ double tolerance,
+ int maxProximalNewtonSteps,
+ double initialRegularization,
+ bool instrumented);
+ void solve();
+
+ void setInitialIterate(const SolutionType& x)
+ {
+ x_ = x;
+ }
+
+ SolutionType getSol() const
+ {
+ return x_;
+ }
+
+ protected:
+#if HAVE_MPI
+ std::unique_ptr<GlobalMapper0> globalMapper0_;
+ std::unique_ptr<GlobalMapper1> globalMapper1_;
+#endif
+ /** \brief The grid */
+ const GridType* grid_;
+
+ /** \brief The solution vectors */
+ SolutionType x_;
+
+ /** \brief The initial regularization parameter for the proximal Newton step */
+ double initialRegularization_;
+ double tolerance_;
+
+ /** \brief Maximum number of proximal-newton steps */
+ std::size_t maxProximalNewtonSteps_;
+
+ /** \brief Hesse matrix */
+ std::unique_ptr<MatrixType> hessianMatrix_;
+
+ /** \brief The assembler for the material law */
+ const MixedGFEAssembler<MixedBasis, TargetSpace0, TargetSpace1>* assembler_;
+
+ /** \brief The solver for the quadratic inner problems */
+ std::shared_ptr<Solvers::CholmodSolver<MatrixType, CorrectionType, BitVector> > innerSolver_;
+ };
+
+} // namespace Dune::GFE
+
+#include "mixedriemannianpnsolver.cc"
+
+#endif
diff --git a/dune/gfe/parallel/mapperfactory.hh b/dune/gfe/parallel/mapperfactory.hh
new file mode 100644
index 0000000000000000000000000000000000000000..4bd2bf23a37651fe543cf90b537e8c4295897bf2
--- /dev/null
+++ b/dune/gfe/parallel/mapperfactory.hh
@@ -0,0 +1,50 @@
+#ifndef DUNE_GFE_PARALLEL_MAPPERFACTORY_HH
+#define DUNE_GFE_PARALLEL_MAPPERFACTORY_HH
+
+#include <dune/grid/utility/globalindexset.hh>
+#include <dune/gfe/parallel/globalmapper.hh>
+#include <dune/gfe/parallel/globalp1mapper.hh>
+#include <dune/gfe/parallel/globalp2mapper.hh>
+#include <dune/gfe/parallel/p2mapper.hh>
+
+namespace Dune::GFE
+{
+
+ /** \brief Assign GlobalMapper and LocalMapper types to a dune-functions function space basis */
+ template <typename Basis>
+ struct MapperFactory
+ {};
+
+ /** \brief Specialization for LagrangeBasis<1> */
+ template <typename GridView>
+ struct MapperFactory<Functions::LagrangeBasis<GridView,1> >
+ {
+ typedef Dune::GlobalP1Mapper<Functions::LagrangeBasis<GridView,1> > GlobalMapper;
+ typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView> LocalMapper;
+ static LocalMapper createLocalMapper(const GridView& gridView)
+ {
+ return LocalMapper(gridView, Dune::mcmgVertexLayout());
+ }
+ };
+
+ template <typename GridView>
+ struct MapperFactory<Dune::Functions::LagrangeBasis<GridView,2> >
+ {
+ typedef Dune::GlobalP2Mapper<Functions::LagrangeBasis<GridView,2> > GlobalMapper;
+ typedef P2BasisMapper<GridView> LocalMapper;
+ static LocalMapper createLocalMapper(const GridView& gridView)
+ {
+ return LocalMapper(gridView);
+ }
+ };
+
+ /** \brief Specialization for LagrangeBasis<3> */
+ template <typename GridView>
+ struct MapperFactory<Dune::Functions::LagrangeBasis<GridView,3> >
+ {
+ // Error: we don't currently have a global P3 mapper
+ };
+
+} // namespace Dune::GFE
+
+#endif
diff --git a/dune/gfe/riemannianpnsolver.cc b/dune/gfe/riemannianpnsolver.cc
index 2637f196b38f4415c6b638f031aeec6b94433d05..836d5b73efd0dbcaeb591a67864a642357fc5aa0 100644
--- a/dune/gfe/riemannianpnsolver.cc
+++ b/dune/gfe/riemannianpnsolver.cc
@@ -101,7 +101,7 @@ setup(const GridType& grid,
operatorAssembler.assembleBulk(Dune::Fufem::istlMatrixBackend(localA), laplaceStiffness);
#if HAVE_MPI
- LocalMapper localMapper = MapperFactory<Basis>::createLocalMapper(grid_->leafGridView());
+ LocalMapper localMapper = Dune::GFE::MapperFactory<Basis>::createLocalMapper(grid_->leafGridView());
MatrixCommunicator<GlobalMapper,
typename GridType::LeafGridView,
@@ -194,7 +194,7 @@ void RiemannianProximalNewtonSolver<Basis,TargetSpace,Assembler>::solve()
VectorCommunicator<GlobalMapper, typename GridType::LeafGridView::Communication, CorrectionType> vectorComm(*globalMapper_,
grid_->leafGridView().comm(),
0);
- LocalMapper localMapper = MapperFactory<Basis>::createLocalMapper(grid_->leafGridView());
+ LocalMapper localMapper = Dune::GFE::MapperFactory<Basis>::createLocalMapper(grid_->leafGridView());
MatrixCommunicator<GlobalMapper,
typename GridType::LeafGridView,
typename GridType::LeafGridView,
diff --git a/dune/gfe/riemannianpnsolver.hh b/dune/gfe/riemannianpnsolver.hh
index 03127e120da75f5167d91bb48ddd564b3edbbbfa..a0a0e7f611a7d115550f1037d7bda0e5c3b76af2 100644
--- a/dune/gfe/riemannianpnsolver.hh
+++ b/dune/gfe/riemannianpnsolver.hh
@@ -16,12 +16,9 @@
#include <dune/solvers/solvers/iterativesolver.hh>
#include <dune/solvers/solvers/cholmodsolver.hh>
-#include "riemanniantrsolver.hh"
-#include <dune/grid/utility/globalindexset.hh>
-#include <dune/gfe/parallel/globalmapper.hh>
-#include <dune/gfe/parallel/globalp1mapper.hh>
-#include <dune/gfe/parallel/globalp2mapper.hh>
-#include <dune/gfe/parallel/p2mapper.hh>
+#include <dune/gfe/assemblers/geodesicfeassembler.hh>
+#include <dune/gfe/parallel/mapperfactory.hh>
+
/** \brief Riemannian proximal-newton solver for geodesic finite-element problems */
template <class Basis, class TargetSpace, class Assembler = GeodesicFEAssembler<Basis,TargetSpace> >
@@ -44,8 +41,8 @@ class RiemannianProximalNewtonSolver
typedef std::vector<TargetSpace> SolutionType;
#if HAVE_MPI
- typedef typename MapperFactory<Basis>::GlobalMapper GlobalMapper;
- typedef typename MapperFactory<Basis>::LocalMapper LocalMapper;
+ typedef typename Dune::GFE::MapperFactory<Basis>::GlobalMapper GlobalMapper;
+ typedef typename Dune::GFE::MapperFactory<Basis>::LocalMapper LocalMapper;
#endif
/** \brief Records information about the last run of the RiemannianProximalNewtonSolver
diff --git a/dune/gfe/riemanniantrsolver.cc b/dune/gfe/riemanniantrsolver.cc
index 941fdd15e3e2830dd648a8826f108ad4dec1dec5..c681055b06ae76e6d621305109c29ce3408bcdb6 100644
--- a/dune/gfe/riemanniantrsolver.cc
+++ b/dune/gfe/riemanniantrsolver.cc
@@ -164,7 +164,7 @@ setup(const GridType& grid,
operatorAssembler.assembleBulk(Dune::Fufem::istlMatrixBackend(localA), laplaceStiffness);
#if HAVE_MPI
- LocalMapper localMapper = MapperFactory<Basis>::createLocalMapper(grid_->leafGridView());
+ LocalMapper localMapper = Dune::GFE::MapperFactory<Basis>::createLocalMapper(grid_->leafGridView());
MatrixCommunicator<GlobalMapper,
typename GridType::LeafGridView,
@@ -384,7 +384,7 @@ void RiemannianTrustRegionSolver<Basis,TargetSpace,Assembler>::solve()
VectorCommunicator<GlobalMapper, typename GridType::LeafGridView::Communication, CorrectionType> vectorComm(*globalMapper_,
grid_->leafGridView().comm(),
0);
- LocalMapper localMapper = MapperFactory<Basis>::createLocalMapper(grid_->leafGridView());
+ LocalMapper localMapper = Dune::GFE::MapperFactory<Basis>::createLocalMapper(grid_->leafGridView());
MatrixCommunicator<GlobalMapper,
typename GridType::LeafGridView,
typename GridType::LeafGridView,
diff --git a/dune/gfe/riemanniantrsolver.hh b/dune/gfe/riemanniantrsolver.hh
index 5893dc357148ace6e59a4f17417f0005bfa52bd7..d2b000ac681a7606701f5dbc2f1cc703e315713a 100644
--- a/dune/gfe/riemanniantrsolver.hh
+++ b/dune/gfe/riemanniantrsolver.hh
@@ -17,46 +17,8 @@
#include <dune/solvers/solvers/loopsolver.hh>
#include <dune/gfe/assemblers/geodesicfeassembler.hh>
-#include <dune/grid/utility/globalindexset.hh>
-#include <dune/gfe/parallel/globalmapper.hh>
-#include <dune/gfe/parallel/globalp1mapper.hh>
-#include <dune/gfe/parallel/globalp2mapper.hh>
-#include <dune/gfe/parallel/p2mapper.hh>
-
-/** \brief Assign GlobalMapper and LocalMapper types to a dune-fufem FunctionSpaceBasis */
-template <typename Basis>
-struct MapperFactory
-{};
-
-/** \brief Specialization for LagrangeBasis<1> */
-template <typename GridView>
-struct MapperFactory<Dune::Functions::LagrangeBasis<GridView,1> >
-{
- typedef Dune::GlobalP1Mapper<Dune::Functions::LagrangeBasis<GridView,1> > GlobalMapper;
- typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView> LocalMapper;
- static LocalMapper createLocalMapper(const GridView& gridView)
- {
- return LocalMapper(gridView, Dune::mcmgVertexLayout());
- }
-};
+#include <dune/gfe/parallel/mapperfactory.hh>
-template <typename GridView>
-struct MapperFactory<Dune::Functions::LagrangeBasis<GridView,2> >
-{
- typedef Dune::GlobalP2Mapper<Dune::Functions::LagrangeBasis<GridView,2> > GlobalMapper;
- typedef P2BasisMapper<GridView> LocalMapper;
- static LocalMapper createLocalMapper(const GridView& gridView)
- {
- return LocalMapper(gridView);
- }
-};
-
-/** \brief Specialization for LagrangeBasis<3> */
-template <typename GridView>
-struct MapperFactory<Dune::Functions::LagrangeBasis<GridView,3> >
-{
- // Error: we don't currently have a global P3 mapper
-};
/** \brief Riemannian trust-region solver for geodesic finite-element problems */
template <class Basis, class TargetSpace, class Assembler = GeodesicFEAssembler<Basis,TargetSpace> >
@@ -79,8 +41,8 @@ class RiemannianTrustRegionSolver
typedef std::vector<TargetSpace> SolutionType;
#if HAVE_MPI
- typedef typename MapperFactory<Basis>::GlobalMapper GlobalMapper;
- typedef typename MapperFactory<Basis>::LocalMapper LocalMapper;
+ typedef typename Dune::GFE::MapperFactory<Basis>::GlobalMapper GlobalMapper;
+ typedef typename Dune::GFE::MapperFactory<Basis>::LocalMapper LocalMapper;
#endif
/** \brief Records information about the last run of the RiemannianTrustRegionSolver
diff --git a/problems/cosserat-continuum-cantilever.parset b/problems/cosserat-continuum-cantilever.parset
index bfa21a8436ebefa2092dc2b059faac5f2d8323fd..ac1a07651c4e668550e5797e222e597df6073367 100644
--- a/problems/cosserat-continuum-cantilever.parset
+++ b/problems/cosserat-continuum-cantilever.parset
@@ -20,6 +20,9 @@ numLevels = 1
# Number of homotopy steps for the Dirichlet boundary conditions
numHomotopySteps = 1
+# Solver type: trustRegion or proximalNewton
+solvertype = trustRegion
+
# Tolerance of the trust region solver
tolerance = 1e-3
diff --git a/problems/film-on-substrate.parset b/problems/film-on-substrate.parset
index 1ce5db4b6676909c3c747f427fddc174dccbc3c7..8bfd20122b0ba089f52c008e4f21f3da2e779291 100644
--- a/problems/film-on-substrate.parset
+++ b/problems/film-on-substrate.parset
@@ -61,7 +61,7 @@ initialDeformation = "x"
# Solver parameters
#############################################
-# Inner solver, cholmod or multigrid
+# Solver type: trustRegion or proximalNewton
solvertype = trustRegion
# Number of homotopy steps for the Dirichlet boundary conditions
diff --git a/src/cosserat-continuum.cc b/src/cosserat-continuum.cc
index ee4f3bd2335f24f577613a8605ecc0cc6762cca1..be9c92887f27d075a9bb6b765eab451d8a1733db 100644
--- a/src/cosserat-continuum.cc
+++ b/src/cosserat-continuum.cc
@@ -43,6 +43,8 @@
#include <dune/foamgrid/foamgrid.hh>
#endif
+#include <dune/istl/multitypeblockvector.hh>
+
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/functions/functionspacebases/compositebasis.hh>
#include <dune/functions/functionspacebases/powerbasis.hh>
@@ -66,6 +68,7 @@
#include <dune/gfe/assemblers/mixedgfeassembler.hh>
#if MIXED_SPACE
+#include <dune/gfe/mixedriemannianpnsolver.hh>
#include <dune/gfe/mixedriemanniantrsolver.hh>
#else
#include <dune/gfe/assemblers/geodesicfeassemblerwrapper.hh>
@@ -504,31 +507,65 @@ int main (int argc, char *argv[]) try
RealTuple<double,3>,
Rotation<double,3> > mixedAssembler(compositeBasis, &localGFEADOLCStiffness);
#if MIXED_SPACE
- MixedRiemannianTrustRegionSolver<GridType,
- CompositeBasis,
- DeformationFEBasis, RealTuple<double,3>,
- OrientationFEBasis, Rotation<double,3> > solver;
- solver.setup(*grid,
- &mixedAssembler,
- deformationFEBasis,
- orientationFEBasis,
- x,
- deformationDirichletDofs,
- orientationDirichletDofs, tolerance,
- maxSolverSteps,
- initialTrustRegionRadius,
- multigridIterations,
- mgTolerance,
- mu, nu1, nu2,
- baseIterations,
- baseTolerance,
- instrumented);
+ if (parameterSet.get<std::string>("solvertype", "trustRegion") == "trustRegion")
+ {
- solver.setScaling(parameterSet.get<FieldVector<double,6> >("solverScaling"));
+ MixedRiemannianTrustRegionSolver<GridType,
+ CompositeBasis,
+ DeformationFEBasis, RealTuple<double,3>,
+ OrientationFEBasis, Rotation<double,3> > solver;
+ solver.setup(*grid,
+ &mixedAssembler,
+ deformationFEBasis,
+ orientationFEBasis,
+ x,
+ deformationDirichletDofs,
+ orientationDirichletDofs, tolerance,
+ maxSolverSteps,
+ initialTrustRegionRadius,
+ multigridIterations,
+ mgTolerance,
+ mu, nu1, nu2,
+ baseIterations,
+ baseTolerance,
+ instrumented);
- solver.setInitialIterate(x);
- solver.solve();
- x = solver.getSol();
+ solver.setScaling(parameterSet.get<FieldVector<double,6> >("solverScaling"));
+
+ solver.setInitialIterate(x);
+ solver.solve();
+ x = solver.getSol();
+ }
+ else
+ {
+ //Create BitVector matching the tangential space
+ const int dimRotationTangent = Rotation<double,3>::TangentVector::dimension;
+ using VectorForBit = MultiTypeBlockVector<std::vector<FieldVector<double,3> >, std::vector<FieldVector<double,dimRotationTangent> > >;
+ using BitVector = Solvers::DefaultBitVector_t<VectorForBit>;
+ BitVector dirichletDofs;
+ dirichletDofs[_0].resize(compositeBasis.size({0}));
+ dirichletDofs[_1].resize(compositeBasis.size({1}));
+ for (size_t i = 0; i < compositeBasis.size({0}); i++) {
+ for (size_t j = 0; j < 3; j++)
+ dirichletDofs[_0][i][j] = deformationDirichletDofs[i][j];
+ }
+ for (size_t i = 0; i < compositeBasis.size({1}); i++) {
+ for (int j = 0; j < dimRotationTangent; j++)
+ dirichletDofs[_1][i][j] = orientationDirichletDofs[i][j];
+ }
+ GFE::MixedRiemannianProximalNewtonSolver<CompositeBasis, DeformationFEBasis, RealTuple<double,3>, OrientationFEBasis, Rotation<double,3>, BitVector> solver;
+ solver.setup(*grid,
+ &mixedAssembler,
+ x,
+ dirichletDofs,
+ tolerance,
+ maxSolverSteps,
+ initialRegularization,
+ instrumented);
+ solver.setInitialIterate(x);
+ solver.solve();
+ x = solver.getSol();
+ }
#else
//The MixedRiemannianTrustRegionSolver can treat the Displacement and Orientation Space as separate ones
//The RiemannianTrustRegionSolver can only treat the Displacement and Rotation together in a ProductManifold.
diff --git a/src/film-on-substrate.cc b/src/film-on-substrate.cc
index 0630bf58bb61cff74976e897171a95bc4248ed0f..97fc8ef8d204ed73bf7a83099e9be19e6980a47b 100644
--- a/src/film-on-substrate.cc
+++ b/src/film-on-substrate.cc
@@ -64,6 +64,7 @@
#include <dune/gfe/assemblers/sumenergy.hh>
#if MIXED_SPACE
+#include <dune/gfe/mixedriemannianpnsolver.hh>
#include <dune/gfe/mixedriemanniantrsolver.hh>
#else
#include <dune/gfe/assemblers/geodesicfeassemblerwrapper.hh>
@@ -650,7 +651,18 @@ int main (int argc, char *argv[]) try
} else { //parameterSet.get<std::string>("solvertype") == "proximalNewton"
#if MIXED_SPACE
- DUNE_THROW(Exception, "Error: There is no MixedRiemannianProximalNewtonSolver!");
+ GFE::MixedRiemannianProximalNewtonSolver<CompositeBasis, DeformationFEBasis, RealTuple<double,dim>, OrientationFEBasis, Rotation<double,dim>, BitVector> solver;
+ solver.setup(*grid,
+ &mixedAssembler,
+ x,
+ dirichletDofs,
+ tolerance,
+ maxSolverSteps,
+ initialRegularization,
+ instrumented);
+ solver.setInitialIterate(x);
+ solver.solve();
+ x = solver.getSol();
#else
RiemannianProximalNewtonSolver<DeformationFEBasis, RBM, GFEAssemblerWrapper> solver;
solver.setup(*grid,
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index ffc3b7de3a65709af4186567953694f9f1c0b114..62f600662aa5d725784d125a0f335efc9b5f48a7 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -53,6 +53,11 @@ dune_add_test(SOURCES geodesicfeassemblerwrappertest.cc
TIMEOUT 600
CMAKE_GUARD MPI_FOUND)
+dune_add_test(SOURCES mixedriemannianpnsolvertest.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/filmonsubstratetest.cc b/test/filmonsubstratetest.cc
index 876871dfe33adb8ae69c98dd06ffce1f65f6cbff..c9b40dda93a72713d1e2cc7b0034606e64856fed 100644
--- a/test/filmonsubstratetest.cc
+++ b/test/filmonsubstratetest.cc
@@ -47,6 +47,7 @@
#include <dune/gfe/assemblers/sumenergy.hh>
#if MIXED_SPACE
+#include <dune/gfe/mixedriemannianpnsolver.hh>
#include <dune/gfe/mixedriemanniantrsolver.hh>
#else
#include <dune/gfe/assemblers/geodesicfeassemblerwrapper.hh>
@@ -458,7 +459,18 @@ int main (int argc, char *argv[])
} else { // solverType == "proximalNewton"
#if MIXED_SPACE
- DUNE_THROW(Exception, "Error: There is no MixedRiemannianProximalNewtonSolver!");
+ GFE::MixedRiemannianProximalNewtonSolver<CompositeBasis, DeformationFEBasis, RealTuple<double,dim>, OrientationFEBasis, Rotation<double,dim>, BitVector> solver;
+ solver.setup(*grid,
+ &mixedAssembler,
+ x,
+ dirichletDofs,
+ tolerance,
+ maxSolverSteps,
+ 1.0 /* initialRegularization */,
+ false /* instrumented */);
+ solver.setInitialIterate(x);
+ solver.solve();
+ x = solver.getSol();
#else
RiemannianProximalNewtonSolver<DeformationFEBasis, RBM, GFEAssemblerWrapper> solver;
solver.setup(*grid,
diff --git a/test/geodesicfeassemblerwrappertest.cc b/test/geodesicfeassemblerwrappertest.cc
index 175cc0b900b8536a2a4f8d3f5687e6ce9b074971..ae8ab8457fafb6adf4c4a58c3afa270f63f43c36 100644
--- a/test/geodesicfeassemblerwrappertest.cc
+++ b/test/geodesicfeassemblerwrappertest.cc
@@ -15,7 +15,6 @@
#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>
@@ -25,9 +24,7 @@
#include <dune/gfe/assemblers/cosseratenergystiffness.hh>
#include <dune/gfe/assemblers/geodesicfeassembler.hh>
#include <dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh>
-#include <dune/gfe/localprojectedfefunction.hh>
#include <dune/gfe/assemblers/mixedgfeassembler.hh>
-#include <dune/gfe/riemanniantrsolver.hh>
#include <dune/gfe/spaces/productmanifold.hh>
#include <dune/gfe/spaces/realtuple.hh>
#include <dune/gfe/spaces/rotation.hh>
@@ -50,8 +47,6 @@ 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> >;
diff --git a/test/mixedriemannianpnsolvertest.cc b/test/mixedriemannianpnsolvertest.cc
new file mode 100644
index 0000000000000000000000000000000000000000..db6b4affb5d126dfa1c82089508a37a78443dee3
--- /dev/null
+++ b/test/mixedriemannianpnsolvertest.cc
@@ -0,0 +1,276 @@
+#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/bitsetvector.hh>
+#include <dune/common/typetraits.hh>
+#include <dune/common/tuplevector.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/fufem/boundarypatch.hh>
+#include <dune/fufem/functiontools/boundarydofs.hh>
+
+#include <dune/grid/utility/structuredgridfactory.hh>
+#include <dune/grid/uggrid.hh>
+
+#include <dune/gfe/assemblers/cosseratenergystiffness.hh>
+#include <dune/gfe/assemblers/geodesicfeassembler.hh>
+#include <dune/gfe/assemblers/geodesicfeassemblerwrapper.hh>
+#include <dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh>
+#include <dune/gfe/assemblers/mixedgfeassembler.hh>
+#include <dune/gfe/mixedriemannianpnsolver.hh>
+#include <dune/gfe/riemannianpnsolver.hh>
+#include <dune/gfe/spaces/productmanifold.hh>
+#include <dune/gfe/spaces/realtuple.hh>
+#include <dune/gfe/spaces/rotation.hh>
+
+/** \file
+ * \brief This test compares the MixedRiemannianPNSolver to the RiemannianPNSolver.
+ */
+
+// grid dimension
+const int gridDim = 2;
+
+// target dimension
+const int dim = 3;
+
+//order of the finite element spaces, they need to be the same to compare!
+const int displacementOrder = 2;
+const int rotationOrder = displacementOrder;
+
+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 = TupleVector<DisplacementVector, RotationVector>;
+using BlockTupleVector = TupleVector<BlockVector<RealTuple<double,dim> >, BlockVector<Rotation<double,dim> > >;
+const int dimRotationTangent = Rotation<double,dim>::TangentVector::dimension;
+
+
+int main (int argc, char *argv[])
+{
+ MPIHelper::instance(argc, argv);
+
+ /////////////////////////////////////////////////////////////////////////
+ // Create the grid
+ /////////////////////////////////////////////////////////////////////////
+ using GridType = UGGrid<gridDim>;
+ auto grid = StructuredGridFactory<GridType>::createCubeGrid({0.0,0.0}, {1.0,1.0}, {2,2});
+ grid->globalRefine(2);
+ grid->loadBalance();
+
+ using GridView = GridType::LeafGridView;
+ GridView gridView = grid->leafGridView();
+
+ /////////////////////////////////////////////////////////////////////////
+ // Create a composite basis and a Lagrange FE basis
+ /////////////////////////////////////////////////////////////////////////
+
+ using namespace Functions::BasisFactory;
+
+ auto compositeBasis = makeBasis(
+ gridView,
+ composite(
+ power<dim>(
+ lagrange<displacementOrder>()
+ ),
+ power<dim>(
+ lagrange<rotationOrder>()
+ )
+ ));
+
+ using CompositeBasis = decltype(compositeBasis);
+
+ using DeformationFEBasis = Functions::LagrangeBasis<GridView,displacementOrder>;
+ DeformationFEBasis deformationFEBasis(gridView);
+
+ /////////////////////////////////////////////////////////////////////////
+ // Create the Neumann and Dirichlet boundary
+ /////////////////////////////////////////////////////////////////////////
+
+ std::function<bool(FieldVector<double,gridDim>)> isNeumann = [](FieldVector<double,gridDim> coordinate) {
+ return coordinate[0] > 0.99; //Neumann for upper boundary
+ };
+ std::function<bool(FieldVector<double,gridDim>)> isDirichlet = [](FieldVector<double,gridDim> coordinate) {
+ return coordinate[0] < 0.01; //Dircichlet for lower boundary
+ };
+
+ BitSetVector<1> neumannVertices(gridView.size(gridDim), false);
+ BitSetVector<1> dirichletVertices(gridView.size(gridDim), false);
+ const GridView::IndexSet& indexSet = gridView.indexSet();
+
+ for (auto&& vertex : vertices(gridView)) {
+ neumannVertices[indexSet.index(vertex)] = isNeumann(vertex.geometry().corner(0));
+ dirichletVertices[indexSet.index(vertex)] = isDirichlet(vertex.geometry().corner(0));
+ }
+
+ BoundaryPatch<GridView> neumannBoundary(gridView, neumannVertices);
+ BoundaryPatch<GridView> dirichletBoundary(gridView, dirichletVertices);
+
+ BitSetVector<1> dirichletNodes(compositeBasis.size({0}),false);
+#if DUNE_VERSION_GTE(DUNE_FUFEM, 2, 10)
+ Fufem::markBoundaryPatchDofs(dirichletBoundary, deformationFEBasis, dirichletNodes);
+#else
+ constructBoundaryDofs(dirichletBoundary, deformationFEBasis, dirichletNodes);
+#endif
+
+ typedef MultiTypeBlockVector<std::vector<FieldVector<double,dim> >, std::vector<FieldVector<double,dimRotationTangent> > > VectorForBit;
+ using BitVector = Solvers::DefaultBitVector_t<VectorForBit>;
+
+ BitVector dirichletDofs;
+ dirichletDofs[_0].resize(compositeBasis.size({0}));
+ dirichletDofs[_1].resize(compositeBasis.size({1}));
+ for (size_t i = 0; i < compositeBasis.size({0}); i++) {
+ for (size_t j = 0; j < dim; j++)
+ dirichletDofs[_0][i][j] = dirichletNodes[i][0];
+ for (size_t j = 0; j < dimRotationTangent; j++)
+ dirichletDofs[_1][i][j] = dirichletNodes[i][0];
+ }
+
+ /////////////////////////////////////////////////////////////////////////
+ // 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";
+
+ parameters["b1"] = "1";
+ parameters["b2"] = "1";
+ parameters["b3"] = "1";
+
+
+ FieldVector<double,dim> values_ = {3e4,2e4,1e4};
+ auto neumannFunction = [&](FieldVector<double, gridDim>){
+ return values_;
+ };
+
+ CosseratEnergyLocalStiffness<decltype(compositeBasis), dim, ValueType> cosseratEnergy(parameters,
+ &neumannBoundary,
+ neumannFunction,
+ nullptr);
+
+ using RBM = GFE::ProductManifold<RealTuple<double,dim>, Rotation<double, dim> >;
+
+ LocalGeodesicFEADOLCStiffness<CompositeBasis,RBM> mixedLocalGFEADOLCStiffness(&cosseratEnergy);
+
+ MixedGFEAssembler<CompositeBasis,
+ RealTuple<double,dim>,
+ Rotation<double,dim> > mixedAssembler(compositeBasis, &mixedLocalGFEADOLCStiffness);
+
+ 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}));
+ BitSetVector<RBM::TangentVector::dimension> dirichletDofsRBM(compositeBasis.size({0}), false);
+
+ for (std::size_t 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];
+ xRBM[i][_0] = x[_0][i];
+ for (int j = 0; j < dim; j ++) { // Displacement part
+ dirichletDofsRBM[i][j] = dirichletDofs[_0][i][j];
+ }
+ xRBM[i][_1] = x[_1][i]; // Rotation part
+ for (int j = dim; j < RBM::TangentVector::dimension; j ++)
+ dirichletDofsRBM[i][j] = dirichletDofs[_1][i][j-dim];
+ }
+
+ //////////////////////////////////////////////////////////////////////////////
+ // Create a MixedRiemannianPNSolver and a normal RiemannianPNSolver
+ // and compare one solver step!
+ //////////////////////////////////////////////////////////////////////////////
+
+ const double tolerance = 1e-7;
+ const int maxSolverSteps = 1;
+ const double initialRegularization = 100;
+ const bool instrumented = false;
+ GFE::MixedRiemannianProximalNewtonSolver<CompositeBasis, DeformationFEBasis, RealTuple<double,dim>, DeformationFEBasis, Rotation<double,dim>, BitVector> mixedSolver;
+ mixedSolver.setup(*grid,
+ &mixedAssembler,
+ x,
+ dirichletDofs,
+ tolerance,
+ maxSolverSteps,
+ initialRegularization,
+ instrumented);
+ mixedSolver.setInitialIterate(x);
+ mixedSolver.solve();
+ x = mixedSolver.getSol();
+
+ RiemannianProximalNewtonSolver<DeformationFEBasis, RBM, GFEAssemblerWrapper> solver;
+ solver.setup(*grid,
+ &assembler,
+ xRBM,
+ dirichletDofsRBM,
+ tolerance,
+ maxSolverSteps,
+ initialRegularization,
+ instrumented);
+ solver.setInitialIterate(xRBM);
+ solver.solve();
+ xRBM = solver.getSol();
+
+ BlockTupleVector xMixed;
+ BlockTupleVector xNotMixed;
+ xNotMixed[_0].resize(compositeBasis.size({0}));
+ xNotMixed[_1].resize(compositeBasis.size({1}));
+ xMixed[_0].resize(compositeBasis.size({0}));
+ xMixed[_1].resize(compositeBasis.size({1}));
+ for (std::size_t i = 0; i < xRBM.size(); i++)
+ {
+ xNotMixed[_0][i] = xRBM[i][_0];
+ xNotMixed[_1][i] = xRBM[i][_1];
+ xMixed[_0][i] = x[_0][i];
+ xMixed[_1][i] = x[_1][i];
+ auto difference0 = xMixed[_0][i].globalCoordinates();
+ auto difference1 = xMixed[_1][i].globalCoordinates();
+ difference0 -= xNotMixed[_0][i].globalCoordinates();
+ difference1 -= xNotMixed[_1][i].globalCoordinates();
+ if (difference0.two_norm() > 1e-1 || difference1.two_norm() > 1e-1) {
+ std::cerr << std::setprecision(9);
+ std::cerr << "At index " << i << " the solution calculated by the MixedRiemannianPNSolver is "
+ << xMixed[_0][i] << " and " << xMixed[_1][i] << " but "
+ << xNotMixed[_0][i] << " and " << xNotMixed[_1][i]
+ << " (calculated by the RiemannianProximalNewtonSolver) was expected!" << std::endl;
+ return 1;
+ }
+ }
+}