Provide access to coarse aggregation for pictures through AMG solver.

Randomize coarsest aggregate numbers to get better contrast between
aggregates.

[[Imported from SVN: r1043]]

istl/paamg/amg.hh

@@ -129,6 +129,13 @@ namespace Dune
       /** \copydoc Preconditioner::post */
       void post(Domain& x);
 
+      /**
+       * @brief Get the aggregate number of each unknown on the coarsest level.
+       * @param cont The random access container to store the numbers in.
+       */
+      template<class A1>
+      void getCoarsestAggregateNumbers(std::vector<std::size_t,A1>& cont);
+
     private:
       /** @brief Multigrid cycle on a level. */
       void mgc(typename Hierarchy<Smoother,A>::Iterator& smoother,
@@ -488,6 +495,14 @@ namespace Dune
       delete &(*rhs_->finest());
       delete rhs_;
     }
+
+    template<class M, class X, class S, class P, class A>
+    template<class A1>
+    void AMG<M,X,S,P,A>::getCoarsestAggregateNumbers(std::vector<std::size_t,A1>& cont)
+    {
+      matrices_->getCoarsestAggregatesOnFinest(cont);
+    }
+
   } // end namespace Amg
 } // end namespace Dune
 

istl/paamg/hierarchy.hh

@@ -397,7 +397,7 @@ namespace Dune
        * @param[out] data The mapping of fine level unknowns to coarse level
        * aggregates.
        */
-      void getCoarsestAggregatesOnFinest(std::vector<std::size_t>& data);
+      void getCoarsestAggregatesOnFinest(std::vector<std::size_t>& data) const;
 
     private:
       typedef typename ConstructionTraits<MatrixOperator>::Arguments MatrixArgs;
@@ -829,7 +829,7 @@ namespace Dune
     }
 
     template<class M, class IS, class A>
-    void MatrixHierarchy<M,IS,A>::getCoarsestAggregatesOnFinest(std::vector<std::size_t>& data)
+    void MatrixHierarchy<M,IS,A>::getCoarsestAggregatesOnFinest(std::vector<std::size_t>& data) const
     {
       int levels=aggregatesMaps().size();
       int maxlevels=parallelInformation_.finest()->communicator().max(levels);
@@ -837,6 +837,7 @@ namespace Dune
       // We need an auxiliary vector for the consecutive prolongation.
       std::vector<std::size_t> tmp;
       std::vector<std::size_t> *coarse, *fine;
+
       // make sure the allocated space suffices.
       tmp.reserve(size);
       data.reserve(size);
@@ -857,14 +858,18 @@ namespace Dune
         std::size_t m=0;
 
         for(typename AggregatesMap::const_iterator iter = map.begin(); iter != map.end(); ++iter)
-          m=std::max(*iter,m);
+          if(*iter< AggregatesMap::ISOLATED)
+            m=std::max(*iter,m);
 
         coarse->resize(m+1);
         std::size_t i=0;
-
+        srand((unsigned)std::clock());
         for(typename std::vector<std::size_t>::iterator iter=coarse->begin(); iter != coarse->end();
             ++iter, ++i)
-          *iter=i;
+        {
+          std::size_t val=(rand());
+          *iter=val;
+        }
       }
 
       typename ParallelInformationHierarchy::Iterator pinfo = parallelInformation().coarsest();