MueLu Version of the Day
Loading...
Searching...
No Matches
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
84namespace MueLu {
85
86 template <class LocalOrdinal, class GlobalOrdinal, class Node>
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: match ML phase2a");
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->getLocalNumElements()
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 */
#define SET_VALID_ENTRY(name)
Container class for aggregation information.
Algorithm for coarsening a graph with uncoupled aggregation.
Among unaggregated points, see if we can make a reasonable size aggregate out of it.
Handle leftover nodes. Try to avoid singleton nodes.
Algorithm for coarsening a graph with structured aggregation.
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
void DeclareInput(Level &currentLevel) const
Input.
void Build(Level &currentLevel) const
Build aggregates.
void BuildInterfaceAggregates(Level &currentLevel, RCP< Aggregates > aggregates, std::vector< unsigned > &aggStat, LO &numNonAggregatedNodes, Array< LO > coarseRate) const
Specifically build aggregates along interfaces.
Algorithm for coarsening a graph with uncoupled aggregation. creates aggregates along an interface us...
Class that holds all level-specific information.
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
int GetLevelID() const
Return level number.
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
static const NoFactory * get()
Algorithm for coarsening a graph with uncoupled aggregation. keep special marked nodes as singleton n...
Builds one-to-one aggregates for all Dirichlet boundary nodes. For some applications this might be ne...
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
Namespace for MueLu classes and methods.
@ Statistics1
Print more statistics.