MueLu  Version of the Day
MueLu_HybridAggregationFactory_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // MueLu: A package for multigrid based preconditioning
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact
39 // Jonathan Hu (jhu@sandia.gov)
40 // Andrey Prokopenko (aprokop@sandia.gov)
41 // Ray Tuminaro (rstumin@sandia.gov)
42 //
43 // ***********************************************************************
44 //
45 // @HEADER
46 #ifndef MUELU_HYBRIDAGGREGATIONFACTORY_DEF_HPP_
47 #define MUELU_HYBRIDAGGREGATIONFACTORY_DEF_HPP_
48 
49 #include <Xpetra_Matrix.hpp>
50 #include <Xpetra_Map.hpp>
51 #include <Xpetra_Vector.hpp>
52 #include <Xpetra_MultiVectorFactory.hpp>
53 #include <Xpetra_VectorFactory.hpp>
54 
56 
57 // Uncoupled Agg
58 #include "MueLu_InterfaceAggregationAlgorithm.hpp"
59 #include "MueLu_OnePtAggregationAlgorithm.hpp"
60 #include "MueLu_PreserveDirichletAggregationAlgorithm.hpp"
61 #include "MueLu_IsolatedNodeAggregationAlgorithm.hpp"
62 
63 #include "MueLu_AggregationPhase1Algorithm.hpp"
64 #include "MueLu_AggregationPhase2aAlgorithm.hpp"
65 #include "MueLu_AggregationPhase2bAlgorithm.hpp"
66 #include "MueLu_AggregationPhase3Algorithm.hpp"
67 
68 // Structured Agg
69 #include "MueLu_AggregationStructuredAlgorithm.hpp"
70 #include "MueLu_UncoupledIndexManager.hpp"
71 //#include "MueLu_LocalLexicographicIndexManager.hpp"
72 //#include "MueLu_GlobalLexicographicIndexManager.hpp"
73 
74 // Shared
75 #include "MueLu_Level.hpp"
76 #include "MueLu_GraphBase.hpp"
77 #include "MueLu_Aggregates.hpp"
78 #include "MueLu_MasterList.hpp"
79 #include "MueLu_Monitor.hpp"
80 #include "MueLu_Utilities.hpp"
81 #include "MueLu_AmalgamationInfo.hpp"
82 
83 
84 namespace MueLu {
85 
86  template <class LocalOrdinal, class GlobalOrdinal, class Node>
88  HybridAggregationFactory() : bDefinitionPhase_(true)
89  { }
90 
91  template <class LocalOrdinal, class GlobalOrdinal, class Node>
94  RCP<ParameterList> validParamList = rcp(new ParameterList());
95 
96  typedef Teuchos::StringToIntegralParameterEntryValidator<int> validatorType;
97 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
98  // From UncoupledAggregationFactory
99  SET_VALID_ENTRY("aggregation: max agg size");
100  SET_VALID_ENTRY("aggregation: min agg size");
101  SET_VALID_ENTRY("aggregation: max selected neighbors");
102  SET_VALID_ENTRY("aggregation: ordering");
103  validParamList->getEntry("aggregation: ordering").setValidator(
104  rcp(new validatorType(Teuchos::tuple<std::string>("natural", "graph", "random"), "aggregation: ordering")));
105  SET_VALID_ENTRY("aggregation: enable phase 1");
106  SET_VALID_ENTRY("aggregation: enable phase 2a");
107  SET_VALID_ENTRY("aggregation: enable phase 2b");
108  SET_VALID_ENTRY("aggregation: enable phase 3");
109  SET_VALID_ENTRY("aggregation: preserve Dirichlet points");
110  SET_VALID_ENTRY("aggregation: allow user-specified singletons");
111  SET_VALID_ENTRY("aggregation: error on nodes with no on-rank neighbors");
112  SET_VALID_ENTRY("aggregation: phase2a include root");
113  SET_VALID_ENTRY("aggregation: phase2a agg factor");
114  SET_VALID_ENTRY("aggregation: phase3 avoid singletons");
115 
116  // From StructuredAggregationFactory
117  SET_VALID_ENTRY("aggregation: coarsening rate");
118  SET_VALID_ENTRY("aggregation: coarsening order");
119  SET_VALID_ENTRY("aggregation: number of spatial dimensions");
120 
121  // From HybridAggregationFactory
122  SET_VALID_ENTRY("aggregation: use interface aggregation");
123 #undef SET_VALID_ENTRY
124 
125  /* From UncoupledAggregation */
126  // general variables needed in AggregationFactory
127  validParamList->set< RCP<const FactoryBase> >("Graph", null, "Generating factory of the graph");
128  validParamList->set< RCP<const FactoryBase> >("DofsPerNode", null, "Generating factory for variable \'DofsPerNode\', usually the same as for \'Graph\'");
129  // special variables necessary for OnePtAggregationAlgorithm
130  validParamList->set<std::string> ("OnePt aggregate map name", "",
131  "Name of input map for single node aggregates. (default='')");
132  validParamList->set<std::string> ("OnePt aggregate map factory", "",
133  "Generating factory of (DOF) map for single node aggregates.");
134 
135  // InterfaceAggregation parameters
136  validParamList->set<std::string> ("Interface aggregate map name", "",
137  "Name of input map for interface aggregates. (default='')");
138  validParamList->set<std::string> ("Interface aggregate map factory", "",
139  "Generating factory of (DOF) map for interface aggregates.");
140  validParamList->set<RCP<const FactoryBase> > ("interfacesDimensions", Teuchos::null,
141  "Describes the dimensions of all the interfaces on this rank.");
142  validParamList->set<RCP<const FactoryBase> > ("nodeOnInterface", Teuchos::null,
143  "List the LIDs of the nodes on any interface.");
144 
145  /* From StructuredAggregation */
146  // general variables needed in AggregationFactory
147  validParamList->set<RCP<const FactoryBase> >("numDimensions", Teuchos::null,
148  "Number of spatial dimension provided by CoordinatesTransferFactory.");
149  validParamList->set<RCP<const FactoryBase> >("lNodesPerDim", Teuchos::null,
150  "Number of nodes per spatial dimmension provided by CoordinatesTransferFactory.");
151 
152 
153  // Hybrid Aggregation Params
154  validParamList->set<RCP<const FactoryBase> > ("aggregationRegionType", Teuchos::null,
155  "Type of aggregation to use on the region (\"structured\" or \"uncoupled\")");
156 
157  return validParamList;
158  }
159 
160  template <class LocalOrdinal, class GlobalOrdinal, class Node>
162  DeclareInput(Level& currentLevel) const {
163  Input(currentLevel, "Graph");
164 
165  ParameterList pL = GetParameterList();
166 
167 
168 
169  /* StructuredAggregation */
170 
171  // Request the local number of nodes per dimensions
172  if(currentLevel.GetLevelID() == 0) {
173  if(currentLevel.IsAvailable("aggregationRegionType", NoFactory::get())) {
174  currentLevel.DeclareInput("aggregationRegionType", NoFactory::get(), this);
175  } else {
176  TEUCHOS_TEST_FOR_EXCEPTION(!currentLevel.IsAvailable("aggregationRegionType",NoFactory::get()),
178  "Aggregation region type was not provided by the user!");
179  }
180  if(currentLevel.IsAvailable("numDimensions", NoFactory::get())) {
181  currentLevel.DeclareInput("numDimensions", NoFactory::get(), this);
182  } else {
183  TEUCHOS_TEST_FOR_EXCEPTION(currentLevel.IsAvailable("numDimensions", NoFactory::get()),
185  "numDimensions was not provided by the user on level0!");
186  }
187  if(currentLevel.IsAvailable("lNodesPerDim", NoFactory::get())) {
188  currentLevel.DeclareInput("lNodesPerDim", NoFactory::get(), this);
189  } else {
190  TEUCHOS_TEST_FOR_EXCEPTION(currentLevel.IsAvailable("lNodesPerDim", NoFactory::get()),
192  "lNodesPerDim was not provided by the user on level0!");
193  }
194  } else {
195  Input(currentLevel, "aggregationRegionType");
196  Input(currentLevel, "numDimensions");
197  Input(currentLevel, "lNodesPerDim");
198  }
199 
200 
201 
202  /* UncoupledAggregation */
203  Input(currentLevel, "DofsPerNode");
204 
205  // request special data necessary for InterfaceAggregation
206  if (pL.get<bool>("aggregation: use interface aggregation") == true){
207  if(currentLevel.GetLevelID() == 0) {
208  if(currentLevel.IsAvailable("interfacesDimensions", NoFactory::get())) {
209  currentLevel.DeclareInput("interfacesDimensions", NoFactory::get(), this);
210  } else {
211  TEUCHOS_TEST_FOR_EXCEPTION(!currentLevel.IsAvailable("interfacesDimensions", NoFactory::get()),
213  "interfacesDimensions was not provided by the user on level0!");
214  }
215  if(currentLevel.IsAvailable("nodeOnInterface", NoFactory::get())) {
216  currentLevel.DeclareInput("nodeOnInterface", NoFactory::get(), this);
217  } else {
218  TEUCHOS_TEST_FOR_EXCEPTION(!currentLevel.IsAvailable("nodeOnInterface", NoFactory::get()),
220  "nodeOnInterface was not provided by the user on level0!");
221  }
222  } else {
223  Input(currentLevel, "interfacesDimensions");
224  Input(currentLevel, "nodeOnInterface");
225  }
226  }
227 
228  // request special data necessary for OnePtAggregationAlgorithm
229  std::string mapOnePtName = pL.get<std::string>("OnePt aggregate map name");
230  if (mapOnePtName.length() > 0) {
231  std::string mapOnePtFactName = pL.get<std::string>("OnePt aggregate map factory");
232  if (mapOnePtFactName == "" || mapOnePtFactName == "NoFactory") {
233  currentLevel.DeclareInput(mapOnePtName, NoFactory::get());
234  } else {
235  RCP<const FactoryBase> mapOnePtFact = GetFactory(mapOnePtFactName);
236  currentLevel.DeclareInput(mapOnePtName, mapOnePtFact.get());
237  }
238  }
239  } // DeclareInput()
240 
241  template <class LocalOrdinal, class GlobalOrdinal, class Node>
243  Build(Level &currentLevel) const {
244  FactoryMonitor m(*this, "Build", currentLevel);
245 
246  RCP<Teuchos::FancyOStream> out;
247  if(const char* dbg = std::getenv("MUELU_HYBRIDAGGREGATION_DEBUG")) {
248  out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
249  out->setShowAllFrontMatter(false).setShowProcRank(true);
250  } else {
251  out = Teuchos::getFancyOStream(rcp(new Teuchos::oblackholestream()));
252  }
253 
254  *out << "Entering hybrid aggregation" << std::endl;
255 
256  ParameterList pL = GetParameterList();
257  bDefinitionPhase_ = false; // definition phase is finished, now all aggregation algorithm information is fixed
258 
259  if (pL.get<int>("aggregation: max agg size") == -1)
260  pL.set("aggregation: max agg size", INT_MAX);
261 
262  // define aggregation algorithms
263  RCP<const FactoryBase> graphFact = GetFactory("Graph");
264 
265  // General problem informations are gathered from data stored in the problem matix.
266  RCP<const GraphBase> graph = Get< RCP<GraphBase> >(currentLevel, "Graph");
267  RCP<const Map> fineMap = graph->GetDomainMap();
268  const int myRank = fineMap->getComm()->getRank();
269  const int numRanks = fineMap->getComm()->getSize();
270 
271  out->setProcRankAndSize(graph->GetImportMap()->getComm()->getRank(),
272  graph->GetImportMap()->getComm()->getSize());
273 
274  // Build aggregates
275  RCP<Aggregates> aggregates = rcp(new Aggregates(*graph));
276  aggregates->setObjectLabel("HB");
277 
278  // construct aggStat information
279  const LO numRows = graph->GetNodeNumVertices();
280  std::vector<unsigned> aggStat(numRows, READY);
281 
282  // Get aggregation type for region
283  std::string regionType;
284  if(currentLevel.GetLevelID() == 0) {
285  // On level 0, data is provided by applications and has no associated factory.
286  regionType = currentLevel.Get<std::string>("aggregationRegionType", NoFactory::get());
287  } else {
288  // On level > 0, data is provided directly by generating factories.
289  regionType = Get< std::string >(currentLevel, "aggregationRegionType");
290  }
291 
292  int numDimensions = 0;
293  if(currentLevel.GetLevelID() == 0) {
294  // On level 0, data is provided by applications and has no associated factory.
295  numDimensions = currentLevel.Get<int>("numDimensions", NoFactory::get());
296  } else {
297  // On level > 0, data is provided directly by generating factories.
298  numDimensions = Get<int>(currentLevel, "numDimensions");
299  }
300 
301  // Get the coarsening rate (potentially used for both structured and uncoupled aggregation if interface)
302  std::string coarseningRate = pL.get<std::string>("aggregation: coarsening rate");
303  Teuchos::Array<LO> coarseRate;
304  try {
305  coarseRate = Teuchos::fromStringToArray<LO>(coarseningRate);
306  } catch(const Teuchos::InvalidArrayStringRepresentation& e) {
307  GetOStream(Errors,-1) << " *** \"aggregation: coarsening rate\" must be a string convertible into an array! *** "
308  << std::endl;
309  throw e;
310  }
311  TEUCHOS_TEST_FOR_EXCEPTION((coarseRate.size() > 1) && (coarseRate.size() < numDimensions),
313  "\"aggregation: coarsening rate\" must have at least as many"
314  " components as the number of spatial dimensions in the problem.");
315 
316  algos_.clear();
317  LO numNonAggregatedNodes = numRows;
318  if (regionType == "structured") {
319  // Add AggregationStructuredAlgorithm
320  algos_.push_back(rcp(new AggregationStructuredAlgorithm(graphFact)));
321 
322  // Since we want to operate on nodes and not dof, we need to modify the rowMap in order to
323  // obtain a nodeMap.
324  const int interpolationOrder = pL.get<int>("aggregation: coarsening order");
325  Array<LO> lFineNodesPerDir(3);
326  if(currentLevel.GetLevelID() == 0) {
327  // On level 0, data is provided by applications and has no associated factory.
328  lFineNodesPerDir = currentLevel.Get<Array<LO> >("lNodesPerDim", NoFactory::get());
329  } else {
330  // On level > 0, data is provided directly by generating factories.
331  lFineNodesPerDir = Get<Array<LO> >(currentLevel, "lNodesPerDim");
332  }
333 
334  // Set lFineNodesPerDir to 1 for directions beyond numDimensions
335  for(int dim = numDimensions; dim < 3; ++dim) {
336  lFineNodesPerDir[dim] = 1;
337  }
338 
339  // Now that we have extracted info from the level, create the IndexManager
340  RCP<MueLu::IndexManager<LO,GO,NO> > geoData;
341  geoData = rcp(new MueLu::UncoupledIndexManager<LO,GO,NO>(fineMap->getComm(),
342  false,
343  numDimensions,
344  interpolationOrder,
345  myRank,
346  numRanks,
347  Array<GO>(3, -1),
348  lFineNodesPerDir,
349  coarseRate, false));
350 
351  TEUCHOS_TEST_FOR_EXCEPTION(fineMap->getNodeNumElements()
352  != static_cast<size_t>(geoData->getNumLocalFineNodes()),
354  "The local number of elements in the graph's map is not equal to "
355  "the number of nodes given by: lNodesPerDim!");
356 
357  aggregates->SetIndexManager(geoData);
358  aggregates->SetNumAggregates(geoData->getNumLocalCoarseNodes());
359 
360  Set(currentLevel, "lCoarseNodesPerDim", geoData->getLocalCoarseNodesPerDir());
361 
362  } // end structured aggregation setup
363 
364  if (regionType == "uncoupled"){
365  // Add unstructred aggregation phases
366  algos_.push_back(rcp(new PreserveDirichletAggregationAlgorithm(graphFact)));
367  if (pL.get<bool>("aggregation: use interface aggregation") == true) algos_.push_back(rcp(new InterfaceAggregationAlgorithm (graphFact)));
368  if (pL.get<bool>("aggregation: allow user-specified singletons") == true) algos_.push_back(rcp(new OnePtAggregationAlgorithm (graphFact)));
369  if (pL.get<bool>("aggregation: enable phase 1" ) == true) algos_.push_back(rcp(new AggregationPhase1Algorithm (graphFact)));
370  if (pL.get<bool>("aggregation: enable phase 2a") == true) algos_.push_back(rcp(new AggregationPhase2aAlgorithm (graphFact)));
371  if (pL.get<bool>("aggregation: enable phase 2b") == true) algos_.push_back(rcp(new AggregationPhase2bAlgorithm (graphFact)));
372  if (pL.get<bool>("aggregation: enable phase 3" ) == true) algos_.push_back(rcp(new AggregationPhase3Algorithm (graphFact)));
373 
374  *out << " Build interface aggregates" << std::endl;
375  // interface
376  if (pL.get<bool>("aggregation: use interface aggregation") == true) {
377  BuildInterfaceAggregates(currentLevel, aggregates, aggStat, numNonAggregatedNodes,
378  coarseRate);
379  }
380 
381  *out << "Treat Dirichlet BC" << std::endl;
382  // Dirichlet boundary
383  ArrayRCP<const bool> dirichletBoundaryMap = graph->GetBoundaryNodeMap();
384  if (dirichletBoundaryMap != Teuchos::null)
385  for (LO i = 0; i < numRows; i++)
386  if (dirichletBoundaryMap[i] == true)
387  aggStat[i] = BOUNDARY;
388 
389  // OnePt aggregation
390  std::string mapOnePtName = pL.get<std::string>("OnePt aggregate map name");
391  RCP<Map> OnePtMap = Teuchos::null;
392  if (mapOnePtName.length()) {
393  std::string mapOnePtFactName = pL.get<std::string>("OnePt aggregate map factory");
394  if (mapOnePtFactName == "" || mapOnePtFactName == "NoFactory") {
395  OnePtMap = currentLevel.Get<RCP<Map> >(mapOnePtName, NoFactory::get());
396  } else {
397  RCP<const FactoryBase> mapOnePtFact = GetFactory(mapOnePtFactName);
398  OnePtMap = currentLevel.Get<RCP<Map> >(mapOnePtName, mapOnePtFact.get());
399  }
400  }
401 
402  LO nDofsPerNode = Get<LO>(currentLevel, "DofsPerNode");
403  GO indexBase = graph->GetDomainMap()->getIndexBase();
404  if (OnePtMap != Teuchos::null) {
405  for (LO i = 0; i < numRows; i++) {
406  // reconstruct global row id (FIXME only works for contiguous maps)
407  GO grid = (graph->GetDomainMap()->getGlobalElement(i)-indexBase) * nDofsPerNode + indexBase;
408  for (LO kr = 0; kr < nDofsPerNode; kr++)
409  if (OnePtMap->isNodeGlobalElement(grid + kr))
410  aggStat[i] = ONEPT;
411  }
412  }
413 
414  // Create a fake lCoarseNodesPerDir for CoordinatesTranferFactory
415  Array<LO> lCoarseNodesPerDir(3,-1);
416  Set(currentLevel, "lCoarseNodesPerDim", lCoarseNodesPerDir);
417  } // end uncoupled aggregation setup
418 
419  aggregates->AggregatesCrossProcessors(false); // No coupled aggregation
420 
421  *out << "Run all the algorithms on the local rank" << std::endl;
422  for (size_t a = 0; a < algos_.size(); a++) {
423  std::string phase = algos_[a]->description();
424  SubFactoryMonitor sfm(*this, "Algo \"" + phase + "\"", currentLevel);
425  *out << regionType <<" | Executing phase " << a << std::endl;
426 
427  int oldRank = algos_[a]->SetProcRankVerbose(this->GetProcRankVerbose());
428  algos_[a]->BuildAggregates(pL, *graph, *aggregates, aggStat, numNonAggregatedNodes);
429  algos_[a]->SetProcRankVerbose(oldRank);
430  *out << regionType <<" | Done Executing phase " << a << std::endl;
431  }
432 
433  *out << "Compute statistics on aggregates" << std::endl;
434  aggregates->ComputeAggregateSizes(true/*forceRecompute*/);
435 
436  Set(currentLevel, "Aggregates", aggregates);
437  Set(currentLevel, "numDimensions", numDimensions);
438  Set(currentLevel, "aggregationRegionTypeCoarse", regionType);
439 
440  GetOStream(Statistics1) << aggregates->description() << std::endl;
441  *out << "HybridAggregation done!" << std::endl;
442  }
443 
444  template <class LocalOrdinal, class GlobalOrdinal, class Node>
446  BuildInterfaceAggregates(Level& currentLevel, RCP<Aggregates> aggregates,
447  std::vector<unsigned>& aggStat, LO& numNonAggregatedNodes,
448  Array<LO> coarseRate) const {
449  FactoryMonitor m(*this, "BuildInterfaceAggregates", currentLevel);
450 
451  RCP<Teuchos::FancyOStream> out;
452  if(const char* dbg = std::getenv("MUELU_HYBRIDAGGREGATION_DEBUG")) {
453  out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
454  out->setShowAllFrontMatter(false).setShowProcRank(true);
455  } else {
456  out = Teuchos::getFancyOStream(rcp(new Teuchos::oblackholestream()));
457  }
458 
459  // Extract and format input data for algo
460  if(coarseRate.size() == 1) {coarseRate.resize(3, coarseRate[0]);}
461  ArrayRCP<LO> vertex2AggId = aggregates->GetVertex2AggId()->getDataNonConst(0);
462  ArrayRCP<LO> procWinner = aggregates->GetProcWinner() ->getDataNonConst(0);
463  Array<LO> interfacesDimensions = Get<Array<LO> >(currentLevel, "interfacesDimensions");
464  Array<LO> nodesOnInterfaces = Get<Array<LO> >(currentLevel, "nodeOnInterface");
465  const int numInterfaces = interfacesDimensions.size() / 3;
466  const int myRank = aggregates->GetMap()->getComm()->getRank();
467 
468  // Create coarse level container to gather data on the fly
469  Array<LO> coarseInterfacesDimensions(interfacesDimensions.size());
470  Array<LO> nodesOnCoarseInterfaces;
471  { // Scoping the temporary variables...
472  LO endRate, totalNumCoarseNodes = 0, numCoarseNodes;
473  for(int interfaceIdx = 0; interfaceIdx < numInterfaces; ++interfaceIdx) {
474  numCoarseNodes = 1;
475  for(int dim = 0; dim < 3; ++dim) {
476  endRate = (interfacesDimensions[3*interfaceIdx + dim] - 1) % coarseRate[dim];
477  if(interfacesDimensions[3*interfaceIdx + dim] == 1) {
478  coarseInterfacesDimensions[3*interfaceIdx + dim] = 1;
479  } else {
480  coarseInterfacesDimensions[3*interfaceIdx + dim]
481  = (interfacesDimensions[3*interfaceIdx+dim]-1) / coarseRate[dim] + 2;
482  if(endRate==0){ coarseInterfacesDimensions[3*interfaceIdx + dim]--;}
483  }
484  numCoarseNodes *= coarseInterfacesDimensions[3*interfaceIdx + dim];
485  }
486  totalNumCoarseNodes += numCoarseNodes;
487  }
488  nodesOnCoarseInterfaces.resize(totalNumCoarseNodes, -1);
489  }
490 
491  Array<LO> endRate(3);
492  LO interfaceOffset = 0, aggregateCount = 0, coarseNodeCount = 0;
493  for(int interfaceIdx = 0; interfaceIdx < numInterfaces; ++interfaceIdx) {
494  ArrayView<LO> fineNodesPerDim = interfacesDimensions(3*interfaceIdx, 3);
495  ArrayView<LO> coarseNodesPerDim = coarseInterfacesDimensions(3*interfaceIdx, 3);
496  LO numInterfaceNodes = 1, numCoarseNodes = 1;
497  for(int dim = 0; dim < 3; ++dim) {
498  numInterfaceNodes *= fineNodesPerDim[dim];
499  numCoarseNodes *= coarseNodesPerDim[dim];
500  endRate[dim] = (fineNodesPerDim[dim]-1) % coarseRate[dim];
501  }
502  ArrayView<LO> interfaceNodes = nodesOnInterfaces(interfaceOffset, numInterfaceNodes);
503 
504  interfaceOffset += numInterfaceNodes;
505 
506  LO rem, rate, fineNodeIdx;
507  Array<LO> nodeIJK(3), coarseIJK(3), rootIJK(3);
508  // First find treat coarse nodes as they generate the aggregate IDs
509  // and they might be repeated on multiple interfaces (think corners and edges).
510  for(LO coarseNodeIdx = 0; coarseNodeIdx < numCoarseNodes; ++coarseNodeIdx) {
511  coarseIJK[2] = coarseNodeIdx / (coarseNodesPerDim[0]*coarseNodesPerDim[1]);
512  rem = coarseNodeIdx % (coarseNodesPerDim[0]*coarseNodesPerDim[1]);
513  coarseIJK[1] = rem / coarseNodesPerDim[0];
514  coarseIJK[0] = rem % coarseNodesPerDim[0];
515 
516  for(LO dim = 0; dim < 3; ++dim) {
517  if(coarseIJK[dim] == coarseNodesPerDim[dim] - 1) {
518  nodeIJK[dim] = fineNodesPerDim[dim] - 1;
519  } else {
520  nodeIJK[dim] = coarseIJK[dim]*coarseRate[dim];
521  }
522  }
523  fineNodeIdx = (nodeIJK[2]*fineNodesPerDim[1] + nodeIJK[1])*fineNodesPerDim[0] + nodeIJK[0];
524 
525  if(aggStat[interfaceNodes[fineNodeIdx]] == READY) {
526  vertex2AggId[interfaceNodes[fineNodeIdx]] = aggregateCount;
527  procWinner[interfaceNodes[fineNodeIdx]] = myRank;
528  aggStat[interfaceNodes[fineNodeIdx]] = AGGREGATED;
529  ++aggregateCount;
530  --numNonAggregatedNodes;
531  }
532  nodesOnCoarseInterfaces[coarseNodeCount] = vertex2AggId[interfaceNodes[fineNodeIdx]];
533  ++coarseNodeCount;
534  }
535 
536  // Now loop over all the node on the interface
537  // skip the coarse nodes as they are already aggregated
538  // and find the appropriate aggregate ID for the fine nodes.
539  for(LO nodeIdx = 0; nodeIdx < numInterfaceNodes; ++nodeIdx) {
540 
541  // If the node is already aggregated skip it!
542  if(aggStat[interfaceNodes[nodeIdx]] == AGGREGATED) {continue;}
543 
544  nodeIJK[2] = nodeIdx / (fineNodesPerDim[0]*fineNodesPerDim[1]);
545  rem = nodeIdx % (fineNodesPerDim[0]*fineNodesPerDim[1]);
546  nodeIJK[1] = rem / fineNodesPerDim[0];
547  nodeIJK[0] = rem % fineNodesPerDim[0];
548 
549  for(int dim = 0; dim < 3; ++dim) {
550  coarseIJK[dim] = nodeIJK[dim] / coarseRate[dim];
551  rem = nodeIJK[dim] % coarseRate[dim];
552  if(nodeIJK[dim] < fineNodesPerDim[dim] - endRate[dim]) {
553  rate = coarseRate[dim];
554  } else {
555  rate = endRate[dim];
556  }
557  if(rem > (rate / 2)) {++coarseIJK[dim];}
558  }
559 
560  for(LO dim = 0; dim < 3; ++dim) {
561  if(coarseIJK[dim] == coarseNodesPerDim[dim] - 1) {
562  nodeIJK[dim] = fineNodesPerDim[dim] - 1;
563  } else {
564  nodeIJK[dim] = coarseIJK[dim]*coarseRate[dim];
565  }
566  }
567  fineNodeIdx = (nodeIJK[2]*fineNodesPerDim[1] + nodeIJK[1])*fineNodesPerDim[0] + nodeIJK[0];
568 
569  vertex2AggId[interfaceNodes[nodeIdx]] = vertex2AggId[interfaceNodes[fineNodeIdx]];
570  procWinner[interfaceNodes[nodeIdx]] = myRank;
571  aggStat[interfaceNodes[nodeIdx]] = AGGREGATED;
572  --numNonAggregatedNodes;
573  } // Loop over interface nodes
574  } // Loop over the interfaces
575 
576  // Update aggregates information before subsequent aggregation algorithms
577  aggregates->SetNumAggregates(aggregateCount);
578 
579  // Set coarse data for next level
580  Set(currentLevel, "coarseInterfacesDimensions", coarseInterfacesDimensions);
581  Set(currentLevel, "nodeOnCoarseInterface", nodesOnCoarseInterfaces);
582 
583  } // BuildInterfaceAggregates()
584 
585 } //namespace MueLu
586 
587 
588 #endif /* MUELU_HYBRIDAGGREGATIONFACTORY_DEF_HPP */
Algorithm for coarsening a graph with uncoupled aggregation. keep special marked nodes as singleton n...
Container class for aggregation information.
Timer to be used in factories. Similar to Monitor but with additional timers.
Algorithm for coarsening a graph with structured aggregation.
Print more statistics.
Namespace for MueLu classes and methods.
static const NoFactory * get()
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:76
Algorithm for coarsening a graph with uncoupled aggregation. creates aggregates along an interface us...
Builds one-to-one aggregates for all Dirichlet boundary nodes. For some applications this might be ne...
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
void DeclareInput(Level &currentLevel) const
Input.
void Build(Level &currentLevel) const
Build aggregates.
Among unaggregated points, see if we can make a reasonable size aggregate out of it.IdeaAmong unaggregated points, see if we can make a reasonable size aggregate out of it. We do this by looking at neighbors and seeing how many are unaggregated and on my processor. Loosely, base the number of new aggregates created on the percentage of unaggregated nodes.
void BuildInterfaceAggregates(Level &currentLevel, RCP< Aggregates > aggregates, std::vector< unsigned > &aggStat, LO &numNonAggregatedNodes, Array< LO > coarseRate) const
Specifically build aggregates along interfaces.
Add leftovers to existing aggregatesIdeaIn phase 2b non-aggregated nodes are added to existing aggreg...
Algorithm for coarsening a graph with uncoupled aggregation.
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need&#39;s value has been saved.
Exception throws to report errors in the internal logical of the program.
Handle leftover nodes. Try to avoid singleton nodesIdeaIn phase 3 we try to stick unaggregated nodes ...
#define SET_VALID_ENTRY(name)
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()