From 9cdb73e452d7a7c9972c46928efc2636ef0485a6 Mon Sep 17 00:00:00 2001
From: Lisa Julia Nebel <lisa_julia.nebel@tu-dresden.de>
Date: Wed, 28 Oct 2020 17:41:47 +0100
Subject: [PATCH] Add a MixedRiemannianPNSolver
The MixedRiemannianPNSolver can assemble when different finite element
spaces are used for deformations and rotations (in, e.g., a Cosserat
model). CHOLMOD is used for the inner solver.
This can be much faster than the MixedRiemannianTRSolver.
---
dune/gfe/mixedriemannianpnsolver.cc | 381 ++++++++++++++++++++++++++++
dune/gfe/mixedriemannianpnsolver.hh | 118 +++++++++
test/CMakeLists.txt | 5 +
test/mixedriemannianpnsolvertest.cc | 276 ++++++++++++++++++++
4 files changed, 780 insertions(+)
create mode 100644 dune/gfe/mixedriemannianpnsolver.cc
create mode 100644 dune/gfe/mixedriemannianpnsolver.hh
create mode 100644 test/mixedriemannianpnsolvertest.cc
diff --git a/dune/gfe/mixedriemannianpnsolver.cc b/dune/gfe/mixedriemannianpnsolver.cc
new file mode 100644
index 00000000..eea8d8cb
--- /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 00000000..a0314bbd
--- /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/test/CMakeLists.txt b/test/CMakeLists.txt
index ffc3b7de..62f60066 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/mixedriemannianpnsolvertest.cc b/test/mixedriemannianpnsolvertest.cc
new file mode 100644
index 00000000..db6b4aff
--- /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;
+ }
+ }
+}
--
GitLab