From 37cc7798d40b5eeb3b703a807d6085f61d6398f1 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Wed, 28 Jan 2015 19:59:45 +0000
Subject: [PATCH] Hand over the FE Basis type as a template parameter
This is much cleaner than handing over the Grid only, and then trying
to guess the basis using the preprocessor. After all, the
RiemannianTrustRegionSolver calls the assembler, so it may know the
function space basis.
[[Imported from SVN: r10035]]
---
dune/gfe/riemanniantrsolver.cc | 20 ++++++++++----------
dune/gfe/riemanniantrsolver.hh | 21 ++++++++-------------
src/cosserat-continuum.cc | 2 +-
src/harmonicmaps.cc | 2 +-
4 files changed, 20 insertions(+), 25 deletions(-)
diff --git a/dune/gfe/riemanniantrsolver.cc b/dune/gfe/riemanniantrsolver.cc
index ffd77f1e..43152330 100644
--- a/dune/gfe/riemanniantrsolver.cc
+++ b/dune/gfe/riemanniantrsolver.cc
@@ -26,10 +26,10 @@
#include <dune/gfe/parallel/matrixcommunicator.hh>
#include <dune/gfe/parallel/vectorcommunicator.hh>
-template <class GridType, class TargetSpace>
-void RiemannianTrustRegionSolver<GridType,TargetSpace>::
+template <class Basis, class TargetSpace>
+void RiemannianTrustRegionSolver<Basis,TargetSpace>::
setup(const GridType& grid,
- const GeodesicFEAssembler<BasisType, TargetSpace>* assembler,
+ const GeodesicFEAssembler<Basis, TargetSpace>* assembler,
const SolutionType& x,
const Dune::BitSetVector<blocksize>& dirichletNodes,
double tolerance,
@@ -112,10 +112,10 @@ setup(const GridType& grid,
// Assemble a Laplace matrix to create a norm that's equivalent to the H1-norm
// //////////////////////////////////////////////////////////////////////////////////////
- BasisType basis(grid.leafGridView());
- OperatorAssembler<BasisType,BasisType> operatorAssembler(basis, basis);
+ Basis basis(grid.leafGridView());
+ OperatorAssembler<Basis,Basis> operatorAssembler(basis, basis);
- LaplaceAssembler<GridType, typename BasisType::LocalFiniteElement, typename BasisType::LocalFiniteElement> laplaceStiffness;
+ LaplaceAssembler<GridType, typename Basis::LocalFiniteElement, typename Basis::LocalFiniteElement> laplaceStiffness;
typedef Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> > ScalarMatrixType;
ScalarMatrixType localA;
@@ -148,7 +148,7 @@ setup(const GridType& grid,
// This will be used to monitor the gradient
// //////////////////////////////////////////////////////////////////////////////////////
- MassAssembler<GridType, typename BasisType::LocalFiniteElement, typename BasisType::LocalFiniteElement> massStiffness;
+ MassAssembler<GridType, typename Basis::LocalFiniteElement, typename Basis::LocalFiniteElement> massStiffness;
ScalarMatrixType localMassMatrix;
operatorAssembler.assemble(massStiffness, localMassMatrix);
@@ -191,7 +191,7 @@ setup(const GridType& grid,
TransferOperatorType pkToP1TransferMatrix;
assembleBasisInterpolationMatrix<TransferOperatorType,
P1NodalBasis<typename GridType::LeafGridView,double>,
- BasisType>(pkToP1TransferMatrix,p1Basis,basis);
+ Basis>(pkToP1TransferMatrix,p1Basis,basis);
// If we are on more than 1 processors, join all local transfer matrices on rank 0,
// and construct a single global transfer operator there.
@@ -261,8 +261,8 @@ setup(const GridType& grid,
}
-template <class GridType, class TargetSpace>
-void RiemannianTrustRegionSolver<GridType,TargetSpace>::solve()
+template <class Basis, class TargetSpace>
+void RiemannianTrustRegionSolver<Basis,TargetSpace>::solve()
{
int argc = 0;
char** argv;
diff --git a/dune/gfe/riemanniantrsolver.hh b/dune/gfe/riemanniantrsolver.hh
index 1880f8a4..118707d6 100644
--- a/dune/gfe/riemanniantrsolver.hh
+++ b/dune/gfe/riemanniantrsolver.hh
@@ -51,11 +51,13 @@ struct MapperFactory<P3NodalBasis<GridView,field_type> >
};
/** \brief Riemannian trust-region solver for geodesic finite-element problems */
-template <class GridType, class TargetSpace>
+template <class Basis, class TargetSpace>
class RiemannianTrustRegionSolver
: public IterativeSolver<std::vector<TargetSpace>,
Dune::BitSetVector<TargetSpace::TangentVector::dimension> >
{
+ typedef typename Basis::GridView::Grid GridType;
+
const static int blocksize = TargetSpace::TangentVector::dimension;
const static int gridDim = GridType::dimension;
@@ -68,15 +70,8 @@ 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
- typedef typename MapperFactory<BasisType>::GlobalMapper GlobalMapper;
- typedef typename MapperFactory<BasisType>::LocalMapper LocalMapper;
+ typedef typename MapperFactory<Basis>::GlobalMapper GlobalMapper;
+ typedef typename MapperFactory<Basis>::LocalMapper LocalMapper;
public:
@@ -90,7 +85,7 @@ public:
/** \brief Set up the solver using a monotone multigrid method as the inner solver */
void setup(const GridType& grid,
- const GeodesicFEAssembler<BasisType, TargetSpace>* assembler,
+ const GeodesicFEAssembler<Basis, TargetSpace>* assembler,
const SolutionType& x,
const Dune::BitSetVector<blocksize>& dirichletNodes,
double tolerance,
@@ -109,7 +104,7 @@ public:
{
scaling_(scaling);
}
-
+
void setIgnoreNodes(const Dune::BitSetVector<blocksize>& ignoreNodes)
{
ignoreNodes_ = &ignoreNodes;
@@ -159,7 +154,7 @@ protected:
std::auto_ptr<MatrixType> hessianMatrix_;
/** \brief The assembler for the material law */
- const GeodesicFEAssembler<BasisType, TargetSpace>* assembler_;
+ const GeodesicFEAssembler<Basis, TargetSpace>* assembler_;
/** \brief The solver for the quadratic inner problems */
std::shared_ptr<Solver> innerSolver_;
diff --git a/src/cosserat-continuum.cc b/src/cosserat-continuum.cc
index 1a4bf2d2..f58b4496 100644
--- a/src/cosserat-continuum.cc
+++ b/src/cosserat-continuum.cc
@@ -271,7 +271,7 @@ int main (int argc, char *argv[]) try
// Create a Riemannian trust-region solver
// /////////////////////////////////////////////////
- RiemannianTrustRegionSolver<GridType,TargetSpace> solver;
+ RiemannianTrustRegionSolver<FEBasis,TargetSpace> solver;
solver.setup(*grid,
&assembler,
x,
diff --git a/src/harmonicmaps.cc b/src/harmonicmaps.cc
index 52bd3b8e..df0abdfc 100644
--- a/src/harmonicmaps.cc
+++ b/src/harmonicmaps.cc
@@ -185,7 +185,7 @@ int main (int argc, char *argv[]) try
// Create a Riemannian trust-region solver
// /////////////////////////////////////////////////
- RiemannianTrustRegionSolver<GridType,TargetSpace> solver;
+ RiemannianTrustRegionSolver<FEBasis,TargetSpace> solver;
solver.setup(*grid,
&assembler,
x,
--
GitLab