143 Level ¤tLevel)
const
145 using Dual2Primal_type = std::map<LocalOrdinal, LocalOrdinal>;
147 const ParameterList &pL = GetParameterList();
149 RCP<const Matrix> A = Get<RCP<Matrix>>(currentLevel,
"A");
152 "Number of dual DOFs per node < 0 (default value). Specify a valid \"number of DOFs per dual node\" in the parameter list for the InterfaceAggregationFactory.");
154 RCP<const Aggregates> primalAggregates = Get<RCP<Aggregates>>(currentLevel,
"Aggregates");
155 ArrayRCP<const LocalOrdinal> primalVertex2AggId = primalAggregates->GetVertex2AggId()->getData(0);
158 RCP<Dual2Primal_type> mapNodesDualToPrimal;
160 mapNodesDualToPrimal = currentLevel.
Get<RCP<Dual2Primal_type>>(
"DualNodeID2PrimalNodeID",
NoFactory::get());
162 mapNodesDualToPrimal = Get<RCP<Dual2Primal_type>>(currentLevel,
"DualNodeID2PrimalNodeID");
164 RCP<const Map> operatorRangeMap = A->getRangeMap();
165 const size_t myRank = operatorRangeMap->getComm()->getRank();
167 LocalOrdinal globalNumDualNodes = operatorRangeMap->getGlobalNumElements() / numDofsPerDualNode;
168 LocalOrdinal localNumDualNodes = operatorRangeMap->getLocalNumElements() / numDofsPerDualNode;
170 TEUCHOS_TEST_FOR_EXCEPTION(localNumDualNodes != Teuchos::as<LocalOrdinal>(mapNodesDualToPrimal->size()),
171 std::runtime_error, prefix <<
" MueLu requires the range map and the DualNodeID2PrimalNodeID map to be compatible.");
173 RCP<const Map> dualNodeMap = Teuchos::null;
174 if (numDofsPerDualNode == 1)
175 dualNodeMap = operatorRangeMap;
179 auto comm = operatorRangeMap->getComm();
180 std::vector<GlobalOrdinal> myDualNodes = {};
182 for (
size_t i = 0; i < operatorRangeMap->getLocalNumElements(); i += numDofsPerDualNode)
183 myDualNodes.push_back((operatorRangeMap->getGlobalElement(i) - indexBase) / numDofsPerDualNode + indexBase);
185 dualNodeMap = MapFactory::Build(operatorRangeMap->lib(), globalNumDualNodes, myDualNodes, indexBase, comm);
187 TEUCHOS_TEST_FOR_EXCEPTION(localNumDualNodes != Teuchos::as<LocalOrdinal>(dualNodeMap->getLocalNumElements()),
188 std::runtime_error, prefix <<
" Local number of dual nodes given by user is incompatible to the dual node map.");
190 RCP<Aggregates> dualAggregates = rcp(
new Aggregates(dualNodeMap));
191 dualAggregates->setObjectLabel(
"InterfaceAggregation");
194 dualAggregates->AggregatesCrossProcessors(primalAggregates->AggregatesCrossProcessors());
196 ArrayRCP<LocalOrdinal> dualVertex2AggId = dualAggregates->GetVertex2AggId()->getDataNonConst(0);
197 ArrayRCP<LocalOrdinal> dualProcWinner = dualAggregates->GetProcWinner()->getDataNonConst(0);
199 RCP<Dual2Primal_type> coarseMapNodesDualToPrimal = rcp(
new Dual2Primal_type());
200 RCP<Dual2Primal_type> coarseMapNodesPrimalToDual = rcp(
new Dual2Primal_type());
209 LocalOrdinal localPrimalNodeID = - Teuchos::ScalarTraits<LocalOrdinal>::one();
210 LocalOrdinal currentPrimalAggId = - Teuchos::ScalarTraits<LocalOrdinal>::one();
211 for (
LocalOrdinal localDualNodeID = 0; localDualNodeID < localNumDualNodes; ++localDualNodeID)
214 localPrimalNodeID = (*mapNodesDualToPrimal)[localDualNodeID];
217 currentPrimalAggId = primalVertex2AggId[localPrimalNodeID];
221 if (coarseMapNodesPrimalToDual->count(currentPrimalAggId) == 0)
224 (*coarseMapNodesPrimalToDual)[currentPrimalAggId] = numLocalDualAggregates;
225 (*coarseMapNodesDualToPrimal)[numLocalDualAggregates] = currentPrimalAggId;
226 ++numLocalDualAggregates;
230 dualVertex2AggId[localDualNodeID] = (*coarseMapNodesPrimalToDual)[currentPrimalAggId];
231 dualProcWinner[localDualNodeID] = myRank;
235 dualAggregates->SetNumAggregates(numLocalDualAggregates);
236 Set(currentLevel,
"Aggregates", dualAggregates);
237 Set(currentLevel,
"CoarseDualNodeID2PrimalNodeID", coarseMapNodesDualToPrimal);
238 GetOStream(
Statistics1) << dualAggregates->description() << std::endl;
243 const std::string& prefix,
Level ¤tLevel)
const
245 const GlobalOrdinal GO_ZERO = Teuchos::ScalarTraits<LocalOrdinal>::zero();
246 const GlobalOrdinal GO_ONE = Teuchos::ScalarTraits<GlobalOrdinal>::one();
253 RCP<const Matrix> A01 = Get<RCP<Matrix>>(currentLevel,
"A");
254 RCP<const Aggregates> primalAggregates = Get<RCP<Aggregates>>(currentLevel,
"Aggregates");
255 ArrayRCP<const LocalOrdinal> primalVertex2AggId = primalAggregates->GetVertex2AggId()->getData(0);
257 auto comm = A01->getRowMap()->getComm();
258 const int myRank = comm->getRank();
260 RCP<const Map> primalInterfaceDofRowMap = Teuchos::null;
263 primalInterfaceDofRowMap = currentLevel.
Get<RCP<const Map>>(
"Primal interface DOF map",
NoFactory::get());
265 primalInterfaceDofRowMap = Get<RCP<const Map>>(currentLevel,
"Primal interface DOF map");
267 TEUCHOS_ASSERT(!primalInterfaceDofRowMap.is_null());
269 if (A01->IsView(
"stridedMaps") && rcp_dynamic_cast<const StridedMap>(A01->getRowMap(
"stridedMaps")) != Teuchos::null) {
270 auto stridedRowMap = rcp_dynamic_cast<const StridedMap>(A01->getRowMap(
"stridedMaps"));
271 auto stridedColMap = rcp_dynamic_cast<const StridedMap>(A01->getColMap(
"stridedMaps"));
272 numDofsPerPrimalNode = Teuchos::as<LocalOrdinal>(stridedRowMap->getFixedBlockSize());
273 numDofsPerDualNode = Teuchos::as<LocalOrdinal>(stridedColMap->getFixedBlockSize());
275 if (numDofsPerPrimalNode != numDofsPerDualNode) {
276 GetOStream(
Warnings) <<
"InterfaceAggregation attempts to work with "
277 << numDofsPerPrimalNode <<
" primal DOFs per node and " << numDofsPerDualNode <<
" dual DOFs per node."
278 <<
"Be careful! Algorithm is not well-tested, if number of primal and dual DOFs per node differ." << std::endl;
283 "InterfaceAggregationFactory could not extract the number of primal DOFs per node from striding information. At least, make sure that StridedMap information has actually been provided.");
285 "InterfaceAggregationFactory could not extract the number of dual DOFs per node from striding information. At least, make sure that StridedMap information has actually been provided.");
306 GlobalOrdinal dualDofOffset = A01->getColMap()->getMinAllGlobalIndex();
310 RCP<const Map> dualDofMap = A01->getDomainMap();
312 dualDofMap->getMaxAllGlobalIndex(), dualBlockDim, dualDofOffset, dualDofMap->getIndexBase());
314 dualDofMap->getMinAllGlobalIndex(), dualBlockDim, dualDofOffset, dualDofMap->getIndexBase());
316 GetOStream(
Runtime1) <<
" Dual DOF map: index base = " << dualDofMap->getIndexBase()
317 <<
", block dim = " << dualBlockDim
318 <<
", gid offset = " << dualDofOffset
321 GetOStream(
Runtime1) <<
" [primal / dual] DOFs per node = [" << numDofsPerPrimalNode
322 <<
"/" << numDofsPerDualNode <<
"]" << std::endl;
325 Array<GlobalOrdinal> dualNodeId2primalNodeId(gMaxDualNodeId - gMinDualNodeId + 1, -GO_ONE);
326 Array<GlobalOrdinal> local_dualNodeId2primalNodeId(gMaxDualNodeId - gMinDualNodeId + 1, -GO_ONE);
329 Array<GlobalOrdinal> dualNodeId2primalAggId(gMaxDualNodeId - gMinDualNodeId + 1, -GO_ONE);
330 Array<GlobalOrdinal> local_dualNodeId2primalAggId(gMaxDualNodeId - gMinDualNodeId + 1, -GO_ONE);
332 Array<GlobalOrdinal> dualDofId2primalDofId(primalInterfaceDofRowMap->getGlobalNumElements(), -GO_ONE);
333 Array<GlobalOrdinal> local_dualDofId2primalDofId(primalInterfaceDofRowMap->getGlobalNumElements(), -GO_ONE);
336 const size_t numMyPrimalInterfaceDOFs = primalInterfaceDofRowMap->getLocalNumElements();
337 for (
size_t r = 0; r < numMyPrimalInterfaceDOFs; r += numDofsPerPrimalNode)
339 GlobalOrdinal gPrimalRowId = primalInterfaceDofRowMap->getGlobalElement(r);
341 if (A01->getRowMap()->isNodeGlobalElement(gPrimalRowId))
343 const LocalOrdinal lPrimalRowId = A01->getRowMap()->getLocalElement(gPrimalRowId);
345 const LocalOrdinal lPrimalNodeId = lPrimalRowId / numDofsPerPrimalNode;
346 const LocalOrdinal primalAggId = primalVertex2AggId[lPrimalNodeId];
348 const GlobalOrdinal gDualDofId = A01->getColMap()->getGlobalElement(r);
352 if (local_dualNodeId2primalNodeId[gDualNodeId - gMinDualNodeId] == -GO_ONE) {
353 local_dualNodeId2primalNodeId[gDualNodeId - gMinDualNodeId] = gPrimalNodeId;
354 local_dualNodeId2primalAggId[gDualNodeId - gMinDualNodeId] = primalAggId;
356 GetOStream(
Warnings) <<
"PROC: " << myRank <<
" gDualNodeId " << gDualNodeId <<
" is already connected to primal nodeId "
357 << local_dualNodeId2primalNodeId[gDualNodeId - gMinDualNodeId]
358 <<
". Ignore new dispNodeId: " << gPrimalNodeId << std::endl;
364 const int dualNodeId2primalNodeIdSize = Teuchos::as<int>(local_dualNodeId2primalNodeId.size());
365 Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, dualNodeId2primalNodeIdSize,
366 &local_dualNodeId2primalNodeId[0], &dualNodeId2primalNodeId[0]);
367 Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, dualNodeId2primalNodeIdSize,
368 &local_dualNodeId2primalAggId[0], &dualNodeId2primalAggId[0]);
373 std::vector<GlobalOrdinal> dualNodes;
374 for (
size_t r = 0; r < A01->getDomainMap()->getLocalNumElements(); r++)
379 GlobalOrdinal gDualDofId = A01->getDomainMap()->getGlobalElement(r);
381 dualNodes.push_back(gDualNodeId);
385 dualNodes.erase(std::unique(dualNodes.begin(), dualNodes.end()), dualNodes.end());
388 Teuchos::RCP<const Map> dualNodeMap = MapFactory::Build(A01->getRowMap()->lib(),
389 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(), dualNodes, A01->getRowMap()->getIndexBase(), comm);
392 Teuchos::RCP<Aggregates> dualAggregates = Teuchos::rcp(
new Aggregates(dualNodeMap));
393 dualAggregates->setObjectLabel(
"UC (dual variables)");
396 Teuchos::ArrayRCP<LocalOrdinal> dualVertex2AggId = dualAggregates->GetVertex2AggId()->getDataNonConst(0);
397 Teuchos::ArrayRCP<LocalOrdinal> dualProcWinner = dualAggregates->GetProcWinner()->getDataNonConst(0);
401 std::map<GlobalOrdinal, LocalOrdinal> primalAggId2localDualAggId;
402 for (
size_t lDualNodeID = 0; lDualNodeID < dualNodeMap->getLocalNumElements(); ++lDualNodeID)
404 const GlobalOrdinal gDualNodeId = dualNodeMap->getGlobalElement(lDualNodeID);
405 const GlobalOrdinal primalAggId = dualNodeId2primalAggId[gDualNodeId - gMinDualNodeId];
406 if (primalAggId2localDualAggId.count(primalAggId) == 0)
407 primalAggId2localDualAggId[primalAggId] = nLocalAggregates++;
408 dualVertex2AggId[lDualNodeID] = primalAggId2localDualAggId[primalAggId];
409 dualProcWinner[lDualNodeID] = myRank;
413 const GlobalOrdinal offset = A01->getColMap()->getMinAllGlobalIndex();
418 RCP<Array<LO>> rowTranslation = rcp(
new Array<LO>());
419 RCP<Array<LO>> colTranslation = rcp(
new Array<LO>());
420 const size_t numMyDualNodes = dualNodeMap->getLocalNumElements();
421 for (
size_t lDualNodeID = 0; lDualNodeID < numMyDualNodes; ++lDualNodeID) {
422 for (
LocalOrdinal dof = 0; dof < numDofsPerDualNode; ++dof) {
423 rowTranslation->push_back(lDualNodeID);
424 colTranslation->push_back(lDualNodeID);
428 TEUCHOS_ASSERT(A01->isFillComplete());
430 RCP<AmalgamationInfo> dualAmalgamationInfo = rcp(
new AmalgamationInfo(rowTranslation, colTranslation,
431 A01->getDomainMap(), A01->getDomainMap(), A01->getDomainMap(),
432 fullblocksize, offset, blockid, nStridedOffset, stridedblocksize));
434 dualAggregates->SetNumAggregates(nLocalAggregates);
435 dualAggregates->AggregatesCrossProcessors(primalAggregates->AggregatesCrossProcessors());
437 if (dualAggregates->AggregatesCrossProcessors())
438 GetOStream(
Runtime1) <<
"Interface aggregates cross processor boundaries." << std::endl;
440 GetOStream(
Runtime1) <<
"Interface aggregates do not cross processor boundaries." << std::endl;
442 currentLevel.
Set(
"Aggregates", dualAggregates,
this);
443 currentLevel.
Set(
"UnAmalgamationInfo", dualAmalgamationInfo,
this);