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"
65template <
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;
87template <
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");
122template <
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 <<
"\".")
141template <
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;
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;
241template <
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;
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);
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Container class for aggregation information.
static const GlobalOrdinal DOFGid2NodeId(GlobalOrdinal gid, LocalOrdinal blockSize, const GlobalOrdinal offset, const GlobalOrdinal indexBase)
Translate global (row/column) id to global amalgamation block id.
minimal container class for storing amalgamation information
Exception throws to report invalid user entry.
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.
void Build(Level ¤tLevel) const override
Build aggregates.
void BuildBasedOnNodeMapping(const std::string &prefix, Level ¤tLevel) const
Build dual aggregates based on a given dual-to-primal node mapping.
void DeclareInput(Level ¤tLevel) const override
Specifies the data that this class needs, and the factories that generate that data.
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.
RCP< const ParameterList > GetValidParameterList() const override
Input.
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)....
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
static const NoFactory * get()
Namespace for MueLu classes and methods.
@ Warnings
Print all warning messages.
@ Statistics1
Print more statistics.
@ Runtime1
Description of what is happening (more verbose)