46 #ifndef MUELU_INTERFACEAGGREGATIONFACTORY_DEF_HPP_ 47 #define MUELU_INTERFACEAGGREGATIONFACTORY_DEF_HPP_ 49 #include <Xpetra_Map.hpp> 50 #include <Xpetra_MapFactory.hpp> 51 #include <Xpetra_StridedMap.hpp> 53 #include "MueLu_Aggregates.hpp" 54 #include "MueLu_AmalgamationFactory.hpp" 55 #include "MueLu_AmalgamationInfo.hpp" 56 #include "MueLu_FactoryManager.hpp" 65 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
68 RCP<ParameterList> validParamList = rcp(
new ParameterList());
70 validParamList->set<RCP<const FactoryBase>>(
"A", Teuchos::null,
"Generating factory of A (matrix block related to dual DOFs)");
71 validParamList->set<RCP<const FactoryBase>>(
"Aggregates", Teuchos::null,
"Generating factory of the Aggregates (for block 0,0)");
73 validParamList->set<std::string>(
"Dual/primal mapping strategy",
"vague",
74 "Strategy to represent mapping between dual and primal quantities [node-based, dof-based]");
76 validParamList->set<RCP<const FactoryBase>>(
"DualNodeID2PrimalNodeID", Teuchos::null,
77 "Generating factory of the DualNodeID2PrimalNodeID map as input data in a Moertel-compatible std::map<LO,LO> to map local IDs of dual nodes to local IDs of primal nodes");
78 validParamList->set<
LocalOrdinal>(
"number of DOFs per dual node", -Teuchos::ScalarTraits<LocalOrdinal>::one(),
79 "Number of DOFs per dual node");
81 validParamList->set<RCP<const FactoryBase>>(
"Primal interface DOF map", Teuchos::null,
82 "Generating factory of the primal DOF row map of slave side of the coupling surface");
84 return validParamList;
87 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
90 Input(currentLevel,
"A");
91 Input(currentLevel,
"Aggregates");
93 const ParameterList &pL = GetParameterList();
95 "Strategy for dual/primal mapping not selected. Please select one of the available strategies.")
96 if (pL.get<std::string>(
"Dual/primal mapping strategy") ==
"node-based")
107 Input(currentLevel,
"DualNodeID2PrimalNodeID");
110 else if (pL.get<std::string>(
"Dual/primal mapping strategy") ==
"dof-based")
115 Input(currentLevel,
"Primal interface DOF map");
122 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
125 const std::string prefix =
"MueLu::InterfaceAggregationFactory::Build: ";
130 const ParameterList &pL = GetParameterList();
131 const std::string parameterName =
"Dual/primal mapping strategy";
132 if (pL.get<std::string>(parameterName) ==
"node-based")
133 BuildBasedOnNodeMapping(prefix, currentLevel);
134 else if (pL.get<std::string>(parameterName) ==
"dof-based")
135 BuildBasedOnPrimalInterfaceDofMap(prefix, currentLevel);
138 "MueLu::InterfaceAggregationFactory::Builld(): Unknown strategy for dual/primal mapping. Set a valid value for the parameter \"" << parameterName <<
"\".")
141 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
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;
159 if (currentLevel.GetLevelID() == 0)
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->getNodeNumElements() / 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->getNodeNumElements(); 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->getNodeNumElements()),
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;
241 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
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;
261 if (currentLevel.GetLevelID() == 0) {
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->getNodeNumElements();
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()->getNodeNumElements(); 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->getNodeNumElements(); ++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->getNodeNumElements();
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);
MueLu::DefaultLocalOrdinal LocalOrdinal
Container class for aggregation information.
Timer to be used in factories. Similar to Monitor but with additional timers.
Namespace for MueLu classes and methods.
static const NoFactory * get()
static const GlobalOrdinal DOFGid2NodeId(GlobalOrdinal gid, LocalOrdinal blockSize, const GlobalOrdinal offset, const GlobalOrdinal indexBase)
Translate global (row/column) id to global amalgamation block id.
int GetLevelID() const
Return level number.
void BuildBasedOnNodeMapping(const std::string &prefix, Level ¤tLevel) const
Build dual aggregates based on a given dual-to-primal node mapping.
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Class that holds all level-specific information.
void Build(Level ¤tLevel) const override
Build aggregates.
void BuildBasedOnPrimalInterfaceDofMap(const std::string &prefix, Level ¤tLevel) const
Build dual aggregates based on a given interface row map of the primal and dual problem.
void DeclareInput(Level ¤tLevel) const override
Specifies the data that this class needs, and the factories that generate that data.
RCP< const ParameterList > GetValidParameterList() const override
Input.
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
Exception throws to report errors in the internal logical of the program.
Print all warning messages.
Description of what is happening (more verbose)
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
minimal container class for storing amalgamation information
Exception throws to report invalid user entry.