From b2573d88955bd42f9300c866b9b89c44b9985f0a Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Mon, 21 Sep 2020 11:09:48 +0200
Subject: [PATCH] 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 | 137 ++++++++++++++++++++
dune/elasticity/common/trustregionsolver.hh | 24 +++-
2 files changed, 159 insertions(+), 2 deletions(-)
diff --git a/dune/elasticity/common/trustregionsolver.cc b/dune/elasticity/common/trustregionsolver.cc
index 82fb9c7..ecae2a3 100644
--- a/dune/elasticity/common/trustregionsolver.cc
+++ b/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()
diff --git a/dune/elasticity/common/trustregionsolver.hh b/dune/elasticity/common/trustregionsolver.hh
index 838f2a1..3947b55 100644
--- a/dune/elasticity/common/trustregionsolver.hh
+++ b/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;
--
GitLab