diff --git a/dune.module b/dune.module
index 6b074074841d86b4e65c247e35d70a750166dab4..e22cda06d446d2df33e6e2422b9c903c266961f2 100644
--- a/dune.module
+++ b/dune.module
@@ -8,4 +8,4 @@ Version: svn
Maintainer: oliver.sander@tu-dresden.de
#depending on
Depends: dune-grid dune-uggrid dune-istl dune-localfunctions dune-functions dune-solvers dune-fufem dune-elasticity
-Suggests: dune-foamgrid
+Suggests: dune-foamgrid dune-parmg
diff --git a/dune/gfe/parallel/globalmapper.hh b/dune/gfe/parallel/globalmapper.hh
new file mode 100644
index 0000000000000000000000000000000000000000..554a70fd67c42d9ea71fabb4e9267738be5bbab1
--- /dev/null
+++ b/dune/gfe/parallel/globalmapper.hh
@@ -0,0 +1,68 @@
+#ifndef DUNE_GFE_PARALLEL_GLOBALMAPPER_HH
+#define DUNE_GFE_PARALLEL_GLOBALMAPPER_HH
+
+#if HAVE_DUNE_PARMG
+#include <dune/parmg/parallel/dofmap.hh>
+#include <dune/parmg/parallel/globaldofindex.hh>
+#include <dune/parmg/parallel/datahandle.hh>
+#endif
+
+namespace Dune {
+
+ template <class Basis>
+ class GlobalMapper
+ {
+ public:
+ using GridView = typename Basis::GridView;
+
+ /** \brief The integer number type used for indices */
+ using Index = typename GridView::IndexSet::IndexType;
+#if HAVE_DUNE_PARMG
+ typedef Dune::ParMG::EntitiesToDofsMap<Basis> DofMap;
+ typedef Dune::ParMG::EntityMasterRank<typename Basis::GridView> Master;
+ typedef std::vector<char> DofMaster;
+
+ GlobalMapper(const Basis& basis)
+ {
+ DofMap dofmap = Dune::ParMG::entitiesToDofsMap(&basis);
+ Master m = Dune::ParMG::entityMasterRank(basis.gridView(), [&](int, int codim) -> bool {
+ return dofmap.codimSet().test(codim);
+ });
+ DofMaster dofmaster = Dune::ParMG::dofMaster(basis, m);
+
+ coarseGlobalDof_ = globalDof(basis,m, dofmap);
+
+ // total number of degrees of freedom
+ size_ = coarseGlobalDof_.size;
+ }
+#else
+ GlobalMapper(const Basis& basis)
+ {
+ // Total number of degrees of freedom
+ size_ = basis.size();
+ }
+#endif
+
+ /** \brief Given a local index, retrieve its index globally unique over all processes. */
+ Index index(const int& localIndex) const {
+#if HAVE_DUNE_PARMG
+ return coarseGlobalDof_.globalDof[localIndex];
+#else
+ return localIndex;
+#endif
+ }
+
+ std::size_t size() const
+ {
+ return size_;
+ }
+
+#if HAVE_DUNE_PARMG
+ ParMG::GlobalDof coarseGlobalDof_;
+#endif
+
+ std::size_t size_;
+ };
+
+}
+#endif // DUNE_GFE_PARALLEL_GLOBALMAPPER_HH
diff --git a/dune/gfe/parallel/globalp1mapper.hh b/dune/gfe/parallel/globalp1mapper.hh
index e7da2815b0f19eb91479592c20c2663d525585ff..3b2fa7302cf16daf401dcaa4137ac724a5db98e9 100644
--- a/dune/gfe/parallel/globalp1mapper.hh
+++ b/dune/gfe/parallel/globalp1mapper.hh
@@ -1,81 +1,59 @@
#ifndef DUNE_GFE_PARALLEL_GLOBALP1MAPPER_HH
#define DUNE_GFE_PARALLEL_GLOBALP1MAPPER_HH
-/** \brief Include standard header files. */
-#include <vector>
-#include <iostream>
-#include <map>
-#include <utility>
-
/** include base class functionality for the communication interface */
-#include <dune/grid/common/datahandleif.hh>
#include <dune/grid/common/mcmgmapper.hh>
-// Include Dune header files
-#include <dune/common/version.hh>
+#include <dune/gfe/parallel/globalmapper.hh>
namespace Dune {
- template <class GridView>
- class GlobalP1Mapper
+ template <class Basis>
+ class GlobalP1Mapper : public GlobalMapper<Basis>
{
- typedef MultipleCodimMultipleGeomTypeMapper<GridView> P1BasisMapper;
-
- public:
+ using GridView = typename Basis::GridView;
+ using P1BasisMapper = MultipleCodimMultipleGeomTypeMapper<GridView>;
+ public:
/** \brief The integer number type used for indices */
- typedef typename GridView::IndexSet::IndexType Index;
-
- typedef std::map<Index,Index> IndexMap;
+ using typename GlobalMapper<Basis>::Index;
- GlobalP1Mapper(const GridView& gridView)
- : p1Mapper_(gridView,mcmgVertexLayout())
+ GlobalP1Mapper(const typename Basis::GridView& gridView)
+ : GlobalMapper<Basis>(gridView),
+ p1Mapper_(gridView,mcmgVertexLayout())
{
+#if !HAVE_DUNE_PARMG
const int dim = GridView::dimension;
GlobalIndexSet<GridView> globalVertexIndexSet(gridView,dim);
// total number of degrees of freedom
- size_ = globalVertexIndexSet.size(dim);
+ this->size_ = globalVertexIndexSet.size(dim);
// Determine
for (auto it = gridView.template begin<0>(); it != gridView.template end<0>(); ++it)
{
// Loop over all vertices
-#if DUNE_VERSION_NEWER(DUNE_GRID,2,4)
for (size_t i=0; i<it->subEntities(dim); i++)
-#else
- for (size_t i=0; i<it->template count<dim>(); i++)
-#endif
{
-#if DUNE_VERSION_NEWER(DUNE_GRID,2,4)
int localIndex = p1Mapper_.subIndex(*it, i, dim);
-#else
- int localIndex = p1Mapper_.map(*it, i, dim);
-#endif
int globalIndex = globalVertexIndexSet.subIndex(*it, i, dim);
localGlobalMap_[localIndex] = globalIndex;
}
-
}
-
- }
-
- /** \brief Given a local index, retrieve its index globally unique over all processes. */
- Index index(const int& localIndex) const {
- return localGlobalMap_.find(localIndex)->second;
+#endif
}
template <class Entity>
Index subIndex(const Entity& entity, uint i, uint codim) const
{
-#if DUNE_VERSION_NEWER(DUNE_GRID,2,4)
int localIndex = p1Mapper_.subIndex(entity, i, codim);
+#if HAVE_DUNE_PARMG
+ return this->index(localIndex);
#else
- int localIndex = p1Mapper_.map(entity, i, codim);
-#endif
return localGlobalMap_.find(localIndex)->second;
+#endif
}
template <class Entity>
@@ -87,18 +65,19 @@ namespace Dune {
return true;
}
- unsigned int size() const
- {
- return size_;
- }
-
P1BasisMapper p1Mapper_;
- IndexMap localGlobalMap_;
+#if !HAVE_DUNE_PARMG
+ typedef std::map<Index,Index> IndexMap;
- size_t size_;
+ /** \brief Given a local index, retrieve its index globally unique over all processes. */
+ Index index(const int& localIndex) const {
+ return localGlobalMap_.find(localIndex)->second;
+ }
+ IndexMap localGlobalMap_;
+#endif
};
}
-#endif /* GLOBALUNIQUEINDEX_HH_ */
+#endif /* DUNE_GFE_PARALLEL_GLOBALP1MAPPER_HH */
diff --git a/dune/gfe/parallel/globalp2mapper.hh b/dune/gfe/parallel/globalp2mapper.hh
index 6bfb74428b26bf89f681210a4c65178175a3cac2..437273b96e4b714dadfa1e9f3eb5258843735b74 100644
--- a/dune/gfe/parallel/globalp2mapper.hh
+++ b/dune/gfe/parallel/globalp2mapper.hh
@@ -1,33 +1,31 @@
#ifndef DUNE_GFE_PARALLEL_GLOBALP2MAPPER_HH
#define DUNE_GFE_PARALLEL_GLOBALP2MAPPER_HH
-/** \brief Include standard header files. */
-#include <vector>
-#include <iostream>
-#include <fstream>
-#include <map>
-#include <utility>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
-/** include base class functionality for the communication interface */
-#include <dune/grid/common/datahandleif.hh>
+#include <dune/gfe/parallel/globalmapper.hh>
-#include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <dune/gfe/parallel/globalmapper.hh>
+
+#include <dune/gfe/parallel/globalmapper.hh>
namespace Dune {
- template <class GridView>
- class GlobalP2Mapper
+ template <class Basis>
+ class GlobalP2Mapper : public Dune::GlobalMapper<Basis>
{
- public:
+ using GridView = typename Basis::GridView;
+ using P2BasisMapper = Functions::LagrangeBasis<GridView,2>;
+ public:
/** \brief The integer number type used for indices */
- typedef typename GridView::IndexSet::IndexType Index;
+ using typename GlobalMapper<Basis>::Index;
- typedef std::map<Index,Index> IndexMap;
-
- GlobalP2Mapper(const GridView& gridView)
- : p2Mapper_(gridView)
+ GlobalP2Mapper(const typename Basis::GridView& gridView) :
+ GlobalMapper<Basis>(gridView),
+ p2Mapper_(gridView)
{
+#if !HAVE_DUNE_PARMG
static_assert(GridView::dimension==2, "Only implemented for two-dimensional grids");
if (gridView.size(GeometryTypes::triangle)>1)
@@ -38,7 +36,7 @@ namespace Dune {
GlobalIndexSet<GridView> globalElementIndex(gridView,0);
// total number of degrees of freedom
- size_ = globalVertexIndex.size(2) + globalEdgeIndex.size(1) + globalElementIndex.size(0);
+ this->size_ = globalVertexIndex.size(2) + globalEdgeIndex.size(1) + globalElementIndex.size(0);
auto localView = p2Mapper_.localView();
@@ -77,21 +75,19 @@ namespace Dune {
localGlobalMap_[localIndex] = globalIndex;
}
-
}
-
- }
-
- /** \brief Given a local index, retrieve its index globally unique over all processes. */
- Index index(const int& localIndex) const {
- return localGlobalMap_.find(localIndex)->second;
+#endif
}
template <class Entity>
Index subIndex(const Entity& entity, uint i, uint codim) const
{
int localIndex = p2Mapper_.map(entity, i, codim);
+#if HAVE_DUNE_PARMG
+ return this->index(localIndex);
+#else
return localGlobalMap_.find(localIndex)->second;
+#endif
}
template <class Entity>
@@ -116,22 +112,24 @@ namespace Dune {
if (not dofFound)
return false;
+#if HAVE_DUNE_PARMG
+ result = this->index(localIndex);
+#else
result = localGlobalMap_.find(localIndex)->second;
+#endif
return true;
}
- unsigned int size() const
- {
- return size_;
+#if !HAVE_DUNE_PARMG
+ /** \brief Given a local index, retrieve its index globally unique over all processes. */
+ Index index(const int& localIndex) const {
+ return localGlobalMap_.find(localIndex)->second;
}
- Functions::LagrangeBasis<GridView,2> p2Mapper_;
-
- IndexMap localGlobalMap_;
-
- size_t size_;
+ std::map<Index,Index> localGlobalMap_;
+#endif
+ P2BasisMapper p2Mapper_;
};
-
}
#endif // DUNE_GFE_PARALLEL_GLOBALP2MAPPER_HH
diff --git a/dune/gfe/riemanniantrsolver.cc b/dune/gfe/riemanniantrsolver.cc
index 8feb344f72223c699bc0ad7ca13b35801efaaf1c..0304ca2c1b60d52267c27e0154db59f1c6e8c669 100644
--- a/dune/gfe/riemanniantrsolver.cc
+++ b/dune/gfe/riemanniantrsolver.cc
@@ -78,12 +78,17 @@ setup(const GridType& grid,
QuadraticIPOptSolver<MatrixType,CorrectionType> baseSolver;
baseSolver.setSolverParameter(baseTolerance, 100, NumProc::QUIET);
#else
+ // First create a Gauss-seidel base solver
+ auto baseSolverStep = std::make_shared<TrustRegionGSStep<MatrixType, CorrectionType>>();
+
// Hack: the two-norm may not scale all that well, but it is fast!
- ::LoopSolver<CorrectionType> baseSolver(TrustRegionGSStep<MatrixType, CorrectionType>{},
- baseIterations,
- baseTolerance,
- TwoNorm<CorrectionType>{},
- Solver::QUIET);
+ auto baseNorm = std::make_shared<TwoNorm<CorrectionType>>();
+
+ auto baseSolver = std::make_shared<::LoopSolver<CorrectionType>>(baseSolverStep,
+ baseIterations,
+ baseTolerance,
+ baseNorm,
+ Solver::QUIET);
#endif
#if HAVE_MPI
// Transfer all Dirichlet data to the master processor
@@ -122,7 +127,7 @@ setup(const GridType& grid,
operatorAssembler.assemble(laplaceStiffness, localA);
#if HAVE_MPI
- LocalMapper localMapper = MapperFactory<typename Basis::GridView,Basis>::createLocalMapper(grid_->leafGridView());
+ LocalMapper localMapper = MapperFactory<Basis>::createLocalMapper(grid_->leafGridView());
MatrixCommunicator<GlobalMapper,
typename GridType::LeafGridView,
@@ -224,7 +229,7 @@ setup(const GridType& grid,
#if HAVE_MPI
// If we are on more than 1 processors, join all local transfer matrices on rank 0,
// and construct a single global transfer operator there.
- typedef Dune::GlobalP1Mapper<typename GridType::LeafGridView> GlobalLeafP1Mapper;
+ typedef Dune::GlobalP1Mapper<Dune::Functions::LagrangeBasis<typename Basis::GridView,1>> GlobalLeafP1Mapper;
GlobalLeafP1Mapper p1Index(grid_->leafGridView());
typedef Dune::MultipleCodimMultipleGeomTypeMapper<typename GridType::LeafGridView> LeafP1LocalMapper;
@@ -251,12 +256,13 @@ setup(const GridType& grid,
for (int i=0; i<numLevels-1; i++) {
// Construct the local multigrid transfer matrix
- TruncatedCompressedMGTransfer<CorrectionType>* newTransferOp = new TruncatedCompressedMGTransfer<CorrectionType>;
+ auto newTransferOp = std::make_unique<TruncatedCompressedMGTransfer<CorrectionType>>();
newTransferOp->setup(*grid_,i,i+1);
#if HAVE_MPI
// If we are on more than 1 processors, join all local transfer matrices on rank 0,
// and construct a single global transfer operator there.
- typedef Dune::GlobalP1Mapper<typename GridType::LevelGridView> GlobalLevelP1Mapper;
+ typedef Dune::Functions::LagrangeBasis<typename GridType::LevelGridView, 1> FEBasis;
+ typedef Dune::GlobalP1Mapper<FEBasis> GlobalLevelP1Mapper;
GlobalLevelP1Mapper fineGUIndex(grid_->levelGridView(i+1));
GlobalLevelP1Mapper coarseGUIndex(grid_->levelGridView(i));
@@ -355,8 +361,7 @@ void RiemannianTrustRegionSolver<Basis,TargetSpace>::solve()
VectorCommunicator<GlobalMapper, typename GridType::LeafGridView::CollectiveCommunication, CorrectionType> vectorComm(*globalMapper_,
grid_->leafGridView().comm(),
0);
-
- LocalMapper localMapper = MapperFactory<typename Basis::GridView,Basis>::createLocalMapper(grid_->leafGridView());
+ LocalMapper localMapper = MapperFactory<Basis>::createLocalMapper(grid_->leafGridView());
MatrixCommunicator<GlobalMapper,
typename GridType::LeafGridView,
typename GridType::LeafGridView,
diff --git a/dune/gfe/riemanniantrsolver.hh b/dune/gfe/riemanniantrsolver.hh
index 7fc9ea26fd2f1ffaca0e7f122194d5a245e2fb28..70c69ffa7a518e134910d8a9a23c439205e7bec1 100644
--- a/dune/gfe/riemanniantrsolver.hh
+++ b/dune/gfe/riemanniantrsolver.hh
@@ -19,20 +19,21 @@
#include "geodesicfeassembler.hh"
#include <dune/grid/utility/globalindexset.hh>
+#include <dune/gfe/parallel/globalmapper.hh>
#include <dune/gfe/parallel/globalp1mapper.hh>
#include <dune/gfe/parallel/globalp2mapper.hh>
#include <dune/gfe/parallel/p2mapper.hh>
/** \brief Assign GlobalMapper and LocalMapper types to a dune-fufem FunctionSpaceBasis */
-template <typename GridView, typename Basis>
+template <typename Basis>
struct MapperFactory
{};
/** \brief Specialization for LagrangeBasis<1> */
template <typename GridView>
-struct MapperFactory<GridView, Dune::Functions::LagrangeBasis<GridView,1> >
+struct MapperFactory<Dune::Functions::LagrangeBasis<GridView,1> >
{
- typedef Dune::GlobalP1Mapper<GridView> GlobalMapper;
+ typedef Dune::GlobalP1Mapper<Dune::Functions::LagrangeBasis<GridView,1>> GlobalMapper;
typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView> LocalMapper;
static LocalMapper createLocalMapper(const GridView& gridView)
{
@@ -43,16 +44,16 @@ struct MapperFactory<GridView, Dune::Functions::LagrangeBasis<GridView,1> >
// This case is not going to actually work, but I need the specialization to make
// the sequential code compile.
template <typename GridView>
-struct MapperFactory<GridView, Dune::Functions::Periodic1DPQ1NodalBasis<GridView> >
+struct MapperFactory<Dune::Functions::Periodic1DPQ1NodalBasis<GridView> >
{
- typedef Dune::GlobalP1Mapper<GridView> GlobalMapper;
+ typedef Dune::GlobalP1Mapper<Dune::Functions::Periodic1DPQ1NodalBasis<GridView>> GlobalMapper;
typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView> LocalMapper;
};
template <typename GridView>
-struct MapperFactory<GridView, Dune::Functions::LagrangeBasis<GridView,2> >
+struct MapperFactory<Dune::Functions::LagrangeBasis<GridView,2> >
{
- typedef Dune::GlobalP2Mapper<GridView> GlobalMapper;
+ typedef Dune::GlobalP2Mapper<Dune::Functions::LagrangeBasis<GridView,2>> GlobalMapper;
typedef P2BasisMapper<GridView> LocalMapper;
static LocalMapper createLocalMapper(const GridView& gridView)
{
@@ -62,7 +63,7 @@ struct MapperFactory<GridView, Dune::Functions::LagrangeBasis<GridView,2> >
/** \brief Specialization for LagrangeBasis<3> */
template <typename GridView>
-struct MapperFactory<GridView, Dune::Functions::LagrangeBasis<GridView,3> >
+struct MapperFactory<Dune::Functions::LagrangeBasis<GridView,3> >
{
// Error: we don't currently have a global P3 mapper
};
@@ -88,8 +89,8 @@ class RiemannianTrustRegionSolver
typedef std::vector<TargetSpace> SolutionType;
#if HAVE_MPI
- typedef typename MapperFactory<typename Basis::GridView,Basis>::GlobalMapper GlobalMapper;
- typedef typename MapperFactory<typename Basis::GridView,Basis>::LocalMapper LocalMapper;
+ typedef typename MapperFactory<Basis>::GlobalMapper GlobalMapper;
+ typedef typename MapperFactory<Basis>::LocalMapper LocalMapper;
#endif
/** \brief Records information about the last run of the RiemannianTrustRegionSolver