Allow to set up a TrustRegionSolver with externally given MMGStep

This is needed as soon as the finite element spaces are not simply
Lagrange spaces anymore, because the current TrustRegionSolver
implementation simply hard-codes particular Lagrange-space
prolongation operators.

dune/elasticity/common/trustregionsolver.cc

@@ -344,6 +344,143 @@ setup(const typename BasisType::GridView::Grid& grid,
 #endif
 }
 
+template <class BasisType, class VectorType>
+void TrustRegionSolver<BasisType,VectorType>::
+setup(const typename BasisType::GridView::Grid& grid,
+      std::shared_ptr<Dune::Elasticity::FEAssembler<BasisType,VectorType> > assembler,
+      const VectorType& x,
+      const Dune::BitSetVector<blocksize>& dirichletNodes,
+      double tolerance,
+      std::size_t maxTrustRegionSteps,
+      double initialTrustRegionRadius,
+      std::shared_ptr<MonotoneMGStep<MatrixType, CorrectionType> > mmgStep,
+      double mgTolerance,
+      std::size_t multigridIterations)
+{
+    grid_                     = &grid;
+    assembler_                = assembler.get();  // FIXME: The class data member should be shared_ptr
+    x_                        = x;
+    this->tolerance_          = tolerance;
+    maxTrustRegionSteps_      = maxTrustRegionSteps;
+    initialTrustRegionRadius_ = initialTrustRegionRadius;
+    innerIterations_          = multigridIterations;
+    innerTolerance_           = mgTolerance;
+    ignoreNodes_              = std::make_shared<Dune::BitSetVector<blocksize>>(dirichletNodes);
+
+    const auto dim = VectorType::value_type::dimension;
+
+#if HAVE_DUNE_PARMG
+    DUNE_THROW(Dune::NotImplemented,
+               "TrustRegionSolver with dune-parmg does not support externally given multigrid step!");
+#endif
+    const auto& basis = assembler_->basis_;
+
+    // //////////////////////////////////////////////////////////////////////////////////////
+    //   Assemble a Laplace matrix to create a norm that's equivalent to the H1-norm
+    // //////////////////////////////////////////////////////////////////////////////////////
+
+
+    using ScalarMatrixType = Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> >;
+    ScalarMatrixType localA;
+
+    Dune::Fufem::DuneFunctionsOperatorAssembler<BasisType,BasisType> operatorAssembler(basis,basis);
+
+    using FiniteElement = std::decay_t<decltype(basis.localView().tree().child(0).finiteElement())>;
+    // construct a Fufem Basis Assembler
+    auto laplaceStiffness = LaplaceAssembler<GridType, FiniteElement, FiniteElement, Dune::ScaledIdentityMatrix<double,dim>>();
+    // transform it to a dune-functions assembler
+    auto localAssembler = [&](const auto& element, auto& localMatrix, auto&& trialLocalView, auto&& ansatzLocalView){
+      // the fufem basis assembler assumes a scalar basis but a blocked element matrix!
+      // we can use ScaledIdentityMatrix here
+      Dune::Matrix<Dune::ScaledIdentityMatrix<double,dim>> localBlockedMatrix( localMatrix.N()/dim, localMatrix.M()/dim );
+      laplaceStiffness.assemble(element,
+                                localBlockedMatrix,
+                                trialLocalView.tree().child(0).finiteElement(),
+                                ansatzLocalView.tree().child(0).finiteElement());
+      // convert back to flat matrix (we need only the diagonal entries in the blocks)
+      localMatrix = 0;
+      for(size_t i=0; i<localBlockedMatrix.N(); i++)
+        for(size_t j=0; j<localBlockedMatrix.M(); j++)
+          Dune::MatrixVector::addToDiagonal(localMatrix[i*dim][j*dim], localBlockedMatrix[i][j].scalar());
+    };
+
+
+    auto matrixBackend = Dune::Fufem::istlMatrixBackend(localA);
+    auto patternBuilder = matrixBackend.patternBuilder();
+
+    operatorAssembler.assembleBulkPattern(patternBuilder);
+    patternBuilder.setupMatrix();
+
+    operatorAssembler.assembleBulkEntries(matrixBackend, localAssembler);
+
+    ScalarMatrixType* A = new ScalarMatrixType(localA);
+
+    h1SemiNorm_ = std::make_shared< H1SemiNorm<CorrectionType> >(*A);
+
+    // //////////////////////////////////////////////////////////////////////////////////////
+    //   Assemble a mass matrix to create a norm that's equivalent to the L2-norm
+    //   This will be used to monitor the gradient
+    // //////////////////////////////////////////////////////////////////////////////////////
+
+    ScalarMatrixType localMassMatrix;
+
+    // construct a Fufem Basis Assembler
+    using FiniteElement = std::decay_t<decltype(basis.localView().tree().child(0).finiteElement())>;
+    auto massStiffness = MassAssembler<GridType, FiniteElement, FiniteElement, Dune::ScaledIdentityMatrix<double,dim>>();
+    // transform it to a dune-functions assembler
+    auto localMassAssembler = [&](const auto& element, auto& localMatrix, auto&& trialLocalView, auto&& ansatzLocalView){
+      // the fufem basis assembler assumes a scalar basis but a blocked element matrix!
+      // we can use ScaledIdentityMatrix here
+      Dune::Matrix<Dune::ScaledIdentityMatrix<double,dim>> localBlockedMatrix( localMatrix.N()/dim, localMatrix.M()/dim );
+      massStiffness.assemble(element,
+                             localBlockedMatrix,
+                             trialLocalView.tree().child(0).finiteElement(),
+                             ansatzLocalView.tree().child(0).finiteElement());
+      // convert back to flat matrix (we need only the diagonal entries in the blocks)
+      localMatrix = 0;
+      for(size_t i=0; i<localBlockedMatrix.N(); i++)
+        for(size_t j=0; j<localBlockedMatrix.M(); j++)
+          Dune::MatrixVector::addToDiagonal(localMatrix[i*dim][j*dim], localBlockedMatrix[i][j].scalar());
+    };
+
+    auto massMatrixBackend = Dune::Fufem::istlMatrixBackend(localMassMatrix);
+    auto massPatternBuilder = massMatrixBackend.patternBuilder();
+
+    operatorAssembler.assembleBulkPattern(massPatternBuilder);
+    massPatternBuilder.setupMatrix();
+
+    operatorAssembler.assembleBulkEntries(massMatrixBackend, localMassAssembler);
+
+    ScalarMatrixType* massMatrix = new ScalarMatrixType(localMassMatrix);
+    l2Norm_ = std::make_shared<H1SemiNorm<CorrectionType> >(*massMatrix);
+
+    // ////////////////////////////////////////////////////////////
+    //    Create Hessian matrix and its occupation structure
+    // ////////////////////////////////////////////////////////////
+
+    hessianMatrix_ = std::make_shared<MatrixType>();
+
+    auto hessianBackend = Dune::Fufem::istlMatrixBackend(*hessianMatrix_);
+    auto hessianPatternBuilder = hessianBackend.patternBuilder();
+    operatorAssembler.assembleBulkPattern(hessianPatternBuilder);
+    hessianPatternBuilder.setupMatrix();
+
+#if !HAVE_DUNE_PARMG
+    innerSolver_ = std::make_shared<LoopSolver<CorrectionType> >(mmgStep,
+                                                                 innerIterations_,
+                                                                 innerTolerance_,
+                                                                 h1SemiNorm_,
+                                                                 Solver::REDUCED);
+
+    ////////////////////////////////////////////////////////////
+    //   Create obstacles
+    ////////////////////////////////////////////////////////////
+
+    hasObstacle_.resize(basis.size(), true);
+    mmgStep->setHasObstacles(hasObstacle_);
+#endif
+}
+
 
 template <class BasisType, class VectorType>
 void TrustRegionSolver<BasisType,VectorType>::solve()

dune/elasticity/common/trustregionsolver.hh

@@ -72,10 +72,15 @@ public:
 
     TrustRegionSolver()
         : IterativeSolver<VectorType, Dune::BitSetVector<blocksize> >(0,100,NumProc::FULL),
-          hessianMatrix_(std::shared_ptr<MatrixType>(NULL)), h1SemiNorm_(NULL)
+          damping_(1.0),
+          hessianMatrix_(std::shared_ptr<MatrixType>(NULL)),
+          h1SemiNorm_(NULL)
     {}
 
-    /** \brief Set up the solver using a monotone multigrid method as the inner solver */
+    /** \brief Set up the solver using a monotone multigrid method as the inner solver
+     *
+     * The monotone multigrid solver is created constructed internally.
+     */
     void setup(const GridType& grid,
                std::conditional_t< // do we have a dune-functions basis? -> choose right assembler type
                   Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>(),
@@ -95,6 +100,21 @@ public:
                double baseTolerance,
                double damping);
 
+    /** \brief Set up the solver using a monotone multigrid method as the inner solver
+     *
+     * The monotone multigrid solver is given from the outside.
+     */
+    void setup(const GridType& grid,
+               std::shared_ptr<Dune::Elasticity::FEAssembler<BasisType,VectorType> > assembler,
+               const SolutionType& x,
+               const Dune::BitSetVector<blocksize>& dirichletNodes,
+               double tolerance,
+               std::size_t maxTrustRegionSteps,
+               double initialTrustRegionRadius,
+               std::shared_ptr<MonotoneMGStep<MatrixType, CorrectionType> > mmgStep,
+               double mgTolerance,
+               std::size_t multigridIterations);
+
     void setIgnoreNodes(const Dune::BitSetVector<blocksize>& ignoreNodes)
     {
         ignoreNodes_ = &ignoreNodes;