From 45039213d1a2d4702e97955bba6f86ceae3b7568 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Sat, 19 Nov 2011 17:08:46 +0000
Subject: [PATCH] add support for third-order Lagrange spaces
[[Imported from SVN: r8234]]
---
dune/gfe/riemanniantrsolver.cc | 34 +++++++++++++++-------------------
dune/gfe/riemanniantrsolver.hh | 22 ++++++++++++----------
2 files changed, 27 insertions(+), 29 deletions(-)
diff --git a/dune/gfe/riemanniantrsolver.cc b/dune/gfe/riemanniantrsolver.cc
index 85f32f7b..bce54849 100644
--- a/dune/gfe/riemanniantrsolver.cc
+++ b/dune/gfe/riemanniantrsolver.cc
@@ -12,8 +12,8 @@
#include <dune/solvers/iterationsteps/trustregiongsstep.hh>
#include <dune/solvers/iterationsteps/mmgstep.hh>
#include <dune/solvers/transferoperators/truncatedcompressedmgtransfer.hh>
-#ifdef HIGHER_ORDER
-#include <dune/solvers/transferoperators/p2top1mgtransfer.hh>
+#if defined THIRD_ORDER || defined SECOND_ORDER
+#include <dune/gfe/pktop1mgtransfer.hh>
#endif
#include <dune/solvers/transferoperators/mandelobsrestrictor.hh>
#include <dune/solvers/solvers/iterativesolver.hh>
@@ -26,11 +26,7 @@
template <class GridType, class TargetSpace>
void RiemannianTrustRegionSolver<GridType,TargetSpace>::
setup(const GridType& grid,
-#ifdef HIGHER_ORDER
- const GeodesicFEAssembler<P2NodalBasis<typename GridType::LeafGridView,double>, TargetSpace>* assembler,
-#else
- const GeodesicFEAssembler<P1NodalBasis<typename GridType::LeafGridView,double>, TargetSpace>* assembler,
-#endif
+ const GeodesicFEAssembler<BasisType, TargetSpace>* assembler,
const SolutionType& x,
const Dune::BitSetVector<blocksize>& dirichletNodes,
double tolerance,
@@ -98,11 +94,11 @@ setup(const GridType& grid,
// //////////////////////////////////////////////////////////////////////////////////////
// Assemble a Laplace matrix to create a norm that's equivalent to the H1-norm
// //////////////////////////////////////////////////////////////////////////////////////
- typedef P1NodalBasis<typename GridType::LeafGridView,double> FEBasis;
- FEBasis basis(grid.leafView());
- OperatorAssembler<FEBasis,FEBasis> operatorAssembler(basis, basis);
+ typedef P1NodalBasis<typename GridType::LeafGridView,double> P1Basis;
+ P1Basis p1Basis(grid.leafView());
+ OperatorAssembler<P1Basis,P1Basis> operatorAssembler(p1Basis, p1Basis);
- LaplaceAssembler<GridType, typename FEBasis::LocalFiniteElement, typename FEBasis::LocalFiniteElement> laplaceStiffness;
+ LaplaceAssembler<GridType, typename P1Basis::LocalFiniteElement, typename P1Basis::LocalFiniteElement> laplaceStiffness;
typedef Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> > ScalarMatrixType;
ScalarMatrixType* A = new ScalarMatrixType;
@@ -138,11 +134,9 @@ setup(const GridType& grid,
// //////////////////////////////////////////////////////////
hasObstacle_.resize(numLevels);
-#ifdef HIGHER_ORDER
- P2NodalBasis<typename GridType::LeafGridView,double> p2Basis(grid_->leafView());
- P1NodalBasis<typename GridType::LeafGridView,double> p1Basis(grid_->leafView());
-
- hasObstacle_.back().resize(p2Basis.size(), true);
+#if defined THIRD_ORDER || defined SECOND_ORDER
+ BasisType basis(grid_->leafView());
+ hasObstacle_.back().resize(basis.size(), true);
for (int i=0; i<hasObstacle_.size()-1; i++)
hasObstacle_[i].resize(grid_->size(i+1, gridDim),true);
@@ -159,10 +153,12 @@ setup(const GridType& grid,
mmgStep->mgTransfer_.resize(numLevels-1);
-#ifdef HIGHER_ORDER
+#if defined THIRD_ORDER || defined SECOND_ORDER
if (numLevels>1) {
- P2toP1MGTransfer<CorrectionType>* topTransferOp = new P2toP1MGTransfer<CorrectionType>;
- topTransferOp->setup(p2Basis,p1Basis);
+ P1NodalBasis<typename GridType::LeafGridView,double> p1Basis(grid_->leafView());
+
+ PKtoP1MGTransfer<CorrectionType>* topTransferOp = new PKtoP1MGTransfer<CorrectionType>;
+ topTransferOp->setup(basis,p1Basis);
mmgStep->mgTransfer_.back() = topTransferOp;
for (int i=0; i<mmgStep->mgTransfer_.size()-1; i++){
diff --git a/dune/gfe/riemanniantrsolver.hh b/dune/gfe/riemanniantrsolver.hh
index 2a1d4642..73ff2e90 100644
--- a/dune/gfe/riemanniantrsolver.hh
+++ b/dune/gfe/riemanniantrsolver.hh
@@ -14,6 +14,8 @@
#include <dune/solvers/solvers/loopsolver.hh>
#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
+#include <dune/fufem/functionspacebases/p2nodalbasis.hh>
+#include <dune/fufem/functionspacebases/p3nodalbasis.hh>
#include "geodesicfeassembler.hh"
@@ -35,6 +37,14 @@ class RiemannianTrustRegionSolver
typedef Dune::BlockVector<Dune::FieldVector<field_type, blocksize> > CorrectionType;
typedef std::vector<TargetSpace> SolutionType;
+#ifdef THIRD_ORDER
+ typedef P3NodalBasis<typename GridType::LeafGridView,double> BasisType;
+#elif defined SECOND_ORDER
+ typedef P2NodalBasis<typename GridType::LeafGridView,double> BasisType;
+#else
+ typedef P1NodalBasis<typename GridType::LeafGridView,double> BasisType;
+#endif
+
public:
RiemannianTrustRegionSolver()
@@ -44,11 +54,7 @@ public:
/** \brief Set up the solver using a monotone multigrid method as the inner solver */
void setup(const GridType& grid,
-#ifdef HIGHER_ORDER
- const GeodesicFEAssembler<P2NodalBasis<typename GridType::LeafGridView,double>, TargetSpace>* assembler,
-#else
- const GeodesicFEAssembler<P1NodalBasis<typename GridType::LeafGridView,double>, TargetSpace>* assembler,
-#endif
+ const GeodesicFEAssembler<BasisType, TargetSpace>* assembler,
const SolutionType& x,
const Dune::BitSetVector<blocksize>& dirichletNodes,
double tolerance,
@@ -103,11 +109,7 @@ protected:
std::auto_ptr<MatrixType> hessianMatrix_;
/** \brief The assembler for the material law */
-#ifdef HIGHER_ORDER
- const GeodesicFEAssembler<P2NodalBasis<typename GridType::LeafGridView,double>, TargetSpace>* assembler_;
-#else
- const GeodesicFEAssembler<P1NodalBasis<typename GridType::LeafGridView,double>, TargetSpace>* assembler_;
-#endif
+ const GeodesicFEAssembler<BasisType, TargetSpace>* assembler_;
/** \brief The solver for the quadratic inner problems */
std::shared_ptr<Solver> innerSolver_;
--
GitLab