71 Aggregates& aggregates, std::vector<unsigned>& aggStat,
72 LO& numNonAggregatedNodes)
const {
73 Monitor m(*
this,
"BuildAggregates");
75 RCP<Teuchos::FancyOStream> out;
76 if(
const char* dbg = std::getenv(
"MUELU_STRUCTUREDALGORITHM_DEBUG")) {
77 out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
78 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
80 out = Teuchos::getFancyOStream(rcp(
new Teuchos::oblackholestream()));
84 const bool coupled = geoData->isAggregationCoupled();
85 const bool singleCoarsePoint = geoData->isSingleCoarsePoint();
86 ArrayRCP<LO> vertex2AggId = aggregates.
GetVertex2AggId()->getDataNonConst(0);
87 ArrayRCP<LO> procWinner = aggregates.
GetProcWinner() ->getDataNonConst(0);
88 Array<LO> ghostedCoarseNodeCoarseLIDs;
89 Array<int> ghostedCoarseNodeCoarsePIDs;
90 Array<GO> ghostedCoarseNodeCoarseGIDs;
92 *out <<
"Extract data for ghosted nodes" << std::endl;
93 geoData->getGhostedNodesData(graph.
GetDomainMap(), ghostedCoarseNodeCoarseLIDs,
94 ghostedCoarseNodeCoarsePIDs, ghostedCoarseNodeCoarseGIDs);
97 Array<LO> ghostedIdx(3), coarseIdx(3);
98 LO ghostedCoarseNodeCoarseLID, aggId;
99 *out <<
"Loop over fine nodes and assign them to an aggregate and a rank" << std::endl;
100 for(LO nodeIdx = 0; nodeIdx < geoData->getNumLocalFineNodes(); ++nodeIdx) {
102 geoData->getFineNodeGhostedTuple(nodeIdx, ghostedIdx[0], ghostedIdx[1], ghostedIdx[2]);
104 for(
int dim = 0; dim < 3; ++dim) {
106 && (geoData->getLocalFineNodesInDir(dim) - 1 < geoData->getCoarseningRate(dim))) {
109 coarseIdx[dim] = ghostedIdx[dim] / geoData->getCoarseningRate(dim);
110 rem = ghostedIdx[dim] % geoData->getCoarseningRate(dim);
111 if(ghostedIdx[dim] - geoData->getOffset(dim)
112 < geoData->getLocalFineNodesInDir(dim) - geoData->getCoarseningEndRate(dim)) {
113 rate = geoData->getCoarseningRate(dim);
115 rate = geoData->getCoarseningEndRate(dim);
117 if(rem > (rate / 2)) {++coarseIdx[dim];}
118 if(coupled && (geoData->getStartGhostedCoarseNode(dim)*geoData->getCoarseningRate(dim)
119 > geoData->getStartIndex(dim))) {--coarseIdx[dim];}
123 geoData->getCoarseNodeGhostedLID(coarseIdx[0], coarseIdx[1], coarseIdx[2],
124 ghostedCoarseNodeCoarseLID);
126 aggId = ghostedCoarseNodeCoarseLIDs[ghostedCoarseNodeCoarseLID];
127 vertex2AggId[nodeIdx] = aggId;
128 procWinner[nodeIdx] = ghostedCoarseNodeCoarsePIDs[ghostedCoarseNodeCoarseLID];
130 --numNonAggregatedNodes;
139 RCP<CrsGraph>& myGraph, RCP<const Map>& coarseCoordinatesFineMap,
140 RCP<const Map>& coarseCoordinatesMap)
const {
141 Monitor m(*
this,
"BuildGraphP");
143 RCP<Teuchos::FancyOStream> out;
144 if(
const char* dbg = std::getenv(
"MUELU_STRUCTUREDALGORITHM_DEBUG")) {
145 out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
146 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
148 out = Teuchos::getFancyOStream(rcp(
new Teuchos::oblackholestream()));
151 const bool coupled = geoData->isAggregationCoupled();
154 int numInterpolationPoints = 0;
155 if(geoData->getInterpolationOrder() == 0) {
156 numInterpolationPoints = 1;
157 }
else if(geoData->getInterpolationOrder() == 1) {
159 numInterpolationPoints = 1 << geoData->getNumDimensions();
161 *out <<
"numInterpolationPoints=" << numInterpolationPoints << std::endl;
163 Array<LO> colIndex((geoData->getNumLocalCoarseNodes() + numInterpolationPoints*
164 (geoData->getNumLocalFineNodes() - geoData->getNumLocalCoarseNodes()))*dofsPerNode);
165 Array<size_t> rowPtr(geoData->getNumLocalFineNodes()*dofsPerNode + 1);
167 ArrayRCP<size_t> nnzOnRow(geoData->getNumLocalFineNodes()*dofsPerNode);
169 *out <<
"Compute prolongatorGraph data" << std::endl;
170 if(geoData->getInterpolationOrder() == 0) {
171 ComputeGraphDataConstant(graph, geoData, dofsPerNode, numInterpolationPoints,
172 nnzOnRow, rowPtr, colIndex);
173 }
else if(geoData->getInterpolationOrder() == 1) {
174 ComputeGraphDataLinear(graph, geoData, dofsPerNode, numInterpolationPoints,
175 nnzOnRow, rowPtr, colIndex);
179 RCP<Map> rowMap = MapFactory::Build(graph.
GetDomainMap(), dofsPerNode);
180 RCP<Map> colMap, domainMap;
181 *out <<
"Compute domain and column maps of the CrsGraph" << std::endl;
183 *out <<
"Extract data for ghosted nodes" << std::endl;
184 Array<LO> ghostedCoarseNodeCoarseLIDs;
185 Array<int> ghostedCoarseNodeCoarsePIDs;
186 Array<GO> ghostedCoarseNodeCoarseGIDs;
187 geoData->getGhostedNodesData(graph.
GetDomainMap(), ghostedCoarseNodeCoarseLIDs,
188 ghostedCoarseNodeCoarsePIDs, ghostedCoarseNodeCoarseGIDs);
193 geoData->getNumGlobalCoarseNodes(),
194 ghostedCoarseNodeCoarseGIDs(),
198 LO coarseNodeIdx = 0;
199 Array<GO> coarseNodeCoarseGIDs, coarseNodeFineGIDs;
200 geoData->getCoarseNodesData(graph.
GetDomainMap(), coarseNodeCoarseGIDs, coarseNodeFineGIDs);
201 for(LO nodeIdx = 0; nodeIdx < ghostedCoarseNodeCoarseGIDs.size(); ++nodeIdx) {
202 if(ghostedCoarseNodeCoarsePIDs[nodeIdx] == colMap->getComm()->getRank()) {
203 coarseNodeCoarseGIDs[coarseNodeIdx] = ghostedCoarseNodeCoarseGIDs[nodeIdx];
207 domainMap = MapFactory::Build(graph.
GetDomainMap()->lib(),
208 geoData->getNumGlobalCoarseNodes(),
209 coarseNodeCoarseGIDs(),
212 coarseCoordinatesMap = MapFactory::Build(graph.
GetDomainMap()->lib(),
213 geoData->getNumGlobalCoarseNodes(),
214 coarseNodeCoarseGIDs(),
217 coarseCoordinatesFineMap = MapFactory::Build(graph.
GetDomainMap()->lib(),
218 geoData->getNumGlobalCoarseNodes(),
219 coarseNodeFineGIDs(),
226 Teuchos::OrdinalTraits<GO>::invalid(),
227 geoData->getNumLocalCoarseNodes()*dofsPerNode,
232 Array<GO> coarseNodeCoarseGIDs(geoData->getNumLocalCoarseNodes());
233 Array<GO> coarseNodeFineGIDs(geoData->getNumLocalCoarseNodes());
234 geoData->getCoarseNodesData(graph.
GetDomainMap(), coarseNodeCoarseGIDs, coarseNodeFineGIDs);
235 coarseCoordinatesMap = MapFactory::Build(graph.
GetDomainMap()->lib(),
236 Teuchos::OrdinalTraits<GO>::invalid(),
237 geoData->getNumLocalCoarseNodes(),
240 coarseCoordinatesFineMap = MapFactory::Build(graph.
GetDomainMap()->lib(),
241 Teuchos::OrdinalTraits<GO>::invalid(),
242 coarseNodeFineGIDs(),
247 *out <<
"Call constructor of CrsGraph" << std::endl;
248 myGraph = CrsGraphFactory::Build(rowMap,
252 *out <<
"Fill CrsGraph" << std::endl;
254 for(LO nodeIdx = 0; nodeIdx < geoData->getNumLocalFineNodes(); ++nodeIdx) {
255 for(LO dof = 0; dof < dofsPerNode; ++dof) {
256 rowIdx = nodeIdx*dofsPerNode + dof;
257 myGraph->insertLocalIndices(rowIdx, colIndex(rowPtr[rowIdx], nnzOnRow[rowIdx]) );
261 *out <<
"Call fillComplete on CrsGraph" << std::endl;
262 myGraph->fillComplete(domainMap, rowMap);
263 *out <<
"Prolongator CrsGraph computed" << std::endl;
271 const LO dofsPerNode,
const int ,
272 ArrayRCP<size_t>& nnzOnRow, Array<size_t>& rowPtr,
273 Array<LO>& colIndex)
const {
275 RCP<Teuchos::FancyOStream> out;
276 if(
const char* dbg = std::getenv(
"MUELU_STRUCTUREDALGORITHM_DEBUG")) {
277 out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
278 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
280 out = Teuchos::getFancyOStream(rcp(
new Teuchos::oblackholestream()));
283 Array<LO> ghostedCoarseNodeCoarseLIDs;
284 Array<int> ghostedCoarseNodeCoarsePIDs;
285 Array<GO> ghostedCoarseNodeCoarseGIDs;
286 geoData->getGhostedNodesData(graph.
GetDomainMap(), ghostedCoarseNodeCoarseLIDs,
287 ghostedCoarseNodeCoarsePIDs, ghostedCoarseNodeCoarseGIDs);
289 LO ghostedCoarseNodeCoarseLID, rem, rate;
290 Array<LO> ghostedIdx(3), coarseIdx(3);
291 for(LO nodeIdx = 0; nodeIdx < geoData->getNumLocalFineNodes(); ++nodeIdx) {
294 geoData->getFineNodeGhostedTuple(nodeIdx, ghostedIdx[0], ghostedIdx[1], ghostedIdx[2]);
296 for(
int dim = 0; dim < 3; ++dim) {
297 if(geoData->isSingleCoarsePoint()
298 && (geoData->getLocalFineNodesInDir(dim) - 1 < geoData->getCoarseningRate(dim))) {
301 coarseIdx[dim] = ghostedIdx[dim] / geoData->getCoarseningRate(dim);
302 rem = ghostedIdx[dim] % geoData->getCoarseningRate(dim);
303 if(ghostedIdx[dim] - geoData->getOffset(dim)
304 < geoData->getLocalFineNodesInDir(dim) - geoData->getCoarseningEndRate(dim)) {
305 rate = geoData->getCoarseningRate(dim);
307 rate = geoData->getCoarseningEndRate(dim);
309 if(rem > (rate / 2)) {++coarseIdx[dim];}
310 if( (geoData->getStartGhostedCoarseNode(dim)*geoData->getCoarseningRate(dim)
311 > geoData->getStartIndex(dim)) && geoData->isAggregationCoupled() ) {
317 geoData->getCoarseNodeGhostedLID(coarseIdx[0], coarseIdx[1], coarseIdx[2],
318 ghostedCoarseNodeCoarseLID);
320 for(LO dof = 0; dof < dofsPerNode; ++dof) {
321 nnzOnRow[nodeIdx*dofsPerNode + dof] = 1;
322 rowPtr[nodeIdx*dofsPerNode + dof + 1] = rowPtr[nodeIdx*dofsPerNode + dof] + 1;
323 colIndex[rowPtr[nodeIdx*dofsPerNode + dof]] =
324 ghostedCoarseNodeCoarseLIDs[ghostedCoarseNodeCoarseLID]*dofsPerNode + dof;
334 const LO dofsPerNode,
const int numInterpolationPoints,
335 ArrayRCP<size_t>& nnzOnRow, Array<size_t>& rowPtr,
336 Array<LO>& colIndex)
const {
338 RCP<Teuchos::FancyOStream> out;
339 if(
const char* dbg = std::getenv(
"MUELU_STRUCTUREDALGORITHM_DEBUG")) {
340 out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
341 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
343 out = Teuchos::getFancyOStream(rcp(
new Teuchos::oblackholestream()));
346 const bool coupled = geoData->isAggregationCoupled();
347 const int numDimensions = geoData->getNumDimensions();
348 Array<LO> ghostedIdx(3,0);
349 Array<LO> coarseIdx(3,0);
350 Array<LO> ijkRem(3,0);
351 const LO coarsePointOffset[8][3] = {{0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {1, 1, 0},
352 {0, 0, 1}, {1, 0, 1}, {0, 1, 1}, {1, 1, 1}};
354 for(LO nodeIdx = 0; nodeIdx < geoData->getNumLocalFineNodes(); ++nodeIdx) {
357 geoData->getFineNodeGhostedTuple(nodeIdx, ghostedIdx[0], ghostedIdx[1], ghostedIdx[2]);
358 for(
int dim=0; dim < numDimensions; dim++){
359 coarseIdx[dim] = ghostedIdx[dim] / geoData->getCoarseningRate(dim);
360 ijkRem[dim] = ghostedIdx[dim] % geoData->getCoarseningRate(dim);
362 if (geoData->getStartGhostedCoarseNode(dim)*geoData->getCoarseningRate(dim)
363 > geoData->getStartIndex(dim)) {
367 if(ghostedIdx[dim] == geoData->getLocalFineNodesInDir(dim) - 1) {
368 coarseIdx[dim] = geoData->getLocalCoarseNodesInDir(dim) - 1;
375 bool allCoarse =
true;
376 Array<bool> isCoarse(numDimensions);
377 for(
int dim = 0; dim < numDimensions; ++dim) {
378 isCoarse[dim] =
false;
380 isCoarse[dim] =
true;
383 if( ghostedIdx[dim]-geoData->getOffset(dim) == geoData->getLocalFineNodesInDir(dim)-1 &&
384 geoData->getMeshEdge(dim*2+1) )
385 isCoarse[dim] =
true;
387 if( ghostedIdx[dim]-geoData->getOffset(dim) == geoData->getLocalFineNodesInDir(dim)-1)
388 isCoarse[dim] =
true;
395 LO rowIdx = 0, colIdx = 0;
397 for(LO dof = 0; dof < dofsPerNode; ++dof) {
398 rowIdx = nodeIdx*dofsPerNode + dof;
399 nnzOnRow[rowIdx] = 1;
400 rowPtr[rowIdx + 1] = rowPtr[rowIdx] + 1;
403 geoData->getCoarseNodeGhostedLID(coarseIdx[0], coarseIdx[1], coarseIdx[2], colIdx);
404 colIndex[rowPtr[rowIdx]] = colIdx*dofsPerNode + dof;
408 for(
int dim = 0; dim < numDimensions; ++dim) {
409 if(coarseIdx[dim] == geoData->getGhostedNodesInDir(dim) - 1)
413 for(LO dof = 0; dof < dofsPerNode; ++dof) {
415 rowIdx = nodeIdx*dofsPerNode + dof;
416 nnzOnRow[rowIdx] = Teuchos::as<size_t>( numInterpolationPoints );
417 rowPtr[rowIdx + 1] = rowPtr[rowIdx] + Teuchos::as<LO>(numInterpolationPoints);
419 for(LO interpIdx = 0; interpIdx < numInterpolationPoints; ++interpIdx) {
420 geoData->getCoarseNodeGhostedLID(coarseIdx[0] + coarsePointOffset[interpIdx][0],
421 coarseIdx[1] + coarsePointOffset[interpIdx][1],
422 coarseIdx[2] + coarsePointOffset[interpIdx][2],
424 colIndex[rowPtr[rowIdx] + interpIdx] = colIdx*dofsPerNode + dof;