10 #ifndef MUELU_INTERFACEAGGREGATIONFACTORY_DEF_HPP_
11 #define MUELU_INTERFACEAGGREGATIONFACTORY_DEF_HPP_
13 #include <Xpetra_Map.hpp>
14 #include <Xpetra_MapFactory.hpp>
15 #include <Xpetra_StridedMap.hpp>
17 #include "MueLu_Aggregates.hpp"
18 #include "MueLu_AmalgamationFactory.hpp"
19 #include "MueLu_AmalgamationInfo.hpp"
27 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
31 validParamList->
set<
RCP<const FactoryBase>>(
"A", Teuchos::null,
"Generating factory of A (matrix block related to dual DOFs)");
32 validParamList->
set<
RCP<const FactoryBase>>(
"Aggregates", Teuchos::null,
"Generating factory of the Aggregates (for block 0,0)");
34 validParamList->
set<std::string>(
"Dual/primal mapping strategy",
"vague",
35 "Strategy to represent mapping between dual and primal quantities [node-based, dof-based]");
38 "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");
40 "Number of DOFs per dual node");
43 "Generating factory of the primal DOF row map of slave side of the coupling surface");
45 return validParamList;
48 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
50 Input(currentLevel,
"A");
51 Input(currentLevel,
"Aggregates");
55 if (pL.
get<std::string>(
"Dual/primal mapping strategy") ==
"node-based") {
62 Input(currentLevel,
"DualNodeID2PrimalNodeID");
64 }
else if (pL.
get<std::string>(
"Dual/primal mapping strategy") ==
"dof-based") {
68 Input(currentLevel,
"Primal interface DOF map");
74 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
76 const std::string prefix =
"MueLu::InterfaceAggregationFactory::Build: ";
82 const std::string parameterName =
"Dual/primal mapping strategy";
83 if (pL.
get<std::string>(parameterName) ==
"node-based")
84 BuildBasedOnNodeMapping(prefix, currentLevel);
85 else if (pL.
get<std::string>(parameterName) ==
"dof-based")
86 BuildBasedOnPrimalInterfaceDofMap(prefix, currentLevel);
89 "MueLu::InterfaceAggregationFactory::Builld(): Unknown strategy for dual/primal mapping. Set a valid value for the parameter \"" << parameterName <<
"\".")
92 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
94 Level ¤tLevel)
const {
95 using Dual2Primal_type = std::map<LocalOrdinal, LocalOrdinal>;
102 "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.");
109 if (currentLevel.GetLevelID() == 0)
112 mapNodesDualToPrimal = Get<RCP<Dual2Primal_type>>(currentLevel,
"DualNodeID2PrimalNodeID");
115 const size_t myRank = operatorRangeMap->getComm()->getRank();
117 GlobalOrdinal dualDofOffset = operatorRangeMap->getMinAllGlobalIndex();
119 LocalOrdinal globalNumDualNodes = operatorRangeMap->getGlobalNumElements() / numDofsPerDualNode;
120 LocalOrdinal localNumDualNodes = operatorRangeMap->getLocalNumElements() / numDofsPerDualNode;
122 TEUCHOS_TEST_FOR_EXCEPTION(localNumDualNodes != Teuchos::as<LocalOrdinal>(mapNodesDualToPrimal->size()),
123 std::runtime_error, prefix <<
" MueLu requires the range map and the DualNodeID2PrimalNodeID map to be compatible.");
126 if (numDofsPerDualNode == 1)
127 dualNodeMap = operatorRangeMap;
130 auto comm = operatorRangeMap->getComm();
131 std::vector<GlobalOrdinal> myDualNodes = {};
133 for (
size_t i = 0; i < operatorRangeMap->getLocalNumElements(); i++) {
134 GlobalOrdinal gDualDofId = operatorRangeMap->getGlobalElement(i);
136 myDualNodes.push_back(gDualNodeId);
140 myDualNodes.erase(std::unique(myDualNodes.begin(), myDualNodes.end()), myDualNodes.end());
142 dualNodeMap = MapFactory::Build(operatorRangeMap->lib(), globalNumDualNodes, myDualNodes, indexBase, comm);
144 TEUCHOS_TEST_FOR_EXCEPTION(localNumDualNodes != Teuchos::as<LocalOrdinal>(dualNodeMap->getLocalNumElements()),
145 std::runtime_error, prefix <<
" Local number of dual nodes given by user is incompatible to the dual node map.");
168 for (
LocalOrdinal localDualNodeID = 0; localDualNodeID < localNumDualNodes; ++localDualNodeID) {
170 localPrimalNodeID = (*mapNodesDualToPrimal)[localDualNodeID];
173 currentPrimalAggId = primalVertex2AggId[localPrimalNodeID];
177 if (coarseMapNodesPrimalToDual->count(currentPrimalAggId) == 0) {
179 (*coarseMapNodesPrimalToDual)[currentPrimalAggId] = numLocalDualAggregates;
180 (*coarseMapNodesDualToPrimal)[numLocalDualAggregates] = currentPrimalAggId;
181 ++numLocalDualAggregates;
185 dualVertex2AggId[localDualNodeID] = (*coarseMapNodesPrimalToDual)[currentPrimalAggId];
186 dualProcWinner[localDualNodeID] = myRank;
195 const size_t numMyDualNodes = dualNodeMap->getLocalNumElements();
196 for (
size_t lDualNodeID = 0; lDualNodeID < numMyDualNodes; ++lDualNodeID) {
197 for (
LocalOrdinal dof = 0; dof < numDofsPerDualNode; ++dof) {
198 rowTranslation->push_back(lDualNodeID);
199 colTranslation->push_back(lDualNodeID);
204 operatorRangeMap, operatorRangeMap, operatorRangeMap,
205 numDofsPerDualNode, dualDofOffset, blockid, nStridedOffset, numDofsPerDualNode));
209 Set(currentLevel,
"Aggregates", dualAggregates);
210 Set(currentLevel,
"CoarseDualNodeID2PrimalNodeID", coarseMapNodesDualToPrimal);
211 Set(currentLevel,
"UnAmalgamationInfo", dualAmalgamationInfo);
215 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
217 const std::string &prefix,
Level ¤tLevel)
const {
230 auto comm = A01->getRowMap()->getComm();
231 const int myRank = comm->getRank();
234 if (currentLevel.GetLevelID() == 0) {
238 primalInterfaceDofRowMap = Get<RCP<const Map>>(currentLevel,
"Primal interface DOF map");
241 if (A01->IsView(
"stridedMaps") && rcp_dynamic_cast<
const StridedMap>(A01->getRowMap(
"stridedMaps")) != Teuchos::null) {
242 auto stridedRowMap = rcp_dynamic_cast<
const StridedMap>(A01->getRowMap(
"stridedMaps"));
243 auto stridedColMap = rcp_dynamic_cast<
const StridedMap>(A01->getColMap(
"stridedMaps"));
244 numDofsPerPrimalNode = Teuchos::as<LocalOrdinal>(stridedRowMap->getFixedBlockSize());
245 numDofsPerDualNode = Teuchos::as<LocalOrdinal>(stridedColMap->getFixedBlockSize());
247 if (numDofsPerPrimalNode != numDofsPerDualNode) {
248 GetOStream(
Warnings) <<
"InterfaceAggregation attempts to work with "
249 << numDofsPerPrimalNode <<
" primal DOFs per node and " << numDofsPerDualNode <<
" dual DOFs per node."
250 <<
"Be careful! Algorithm is not well-tested, if number of primal and dual DOFs per node differ." << std::endl;
255 "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.");
257 "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.");
278 GlobalOrdinal dualDofOffset = A01->getRowMap()->getMaxAllGlobalIndex() + 1;
283 dualDofMap->getMaxAllGlobalIndex(), dualBlockDim, dualDofOffset, dualDofMap->getIndexBase());
285 dualDofMap->getMinAllGlobalIndex(), dualBlockDim, dualDofOffset, dualDofMap->getIndexBase());
287 GetOStream(
Runtime1) <<
" Dual DOF map: index base = " << dualDofMap->getIndexBase()
288 <<
", block dim = " << dualBlockDim
289 <<
", gid offset = " << dualDofOffset
292 GetOStream(
Runtime1) <<
" [primal / dual] DOFs per node = [" << numDofsPerPrimalNode
293 <<
"/" << numDofsPerDualNode <<
"]" << std::endl;
297 Array<GlobalOrdinal> local_dualNodeId2primalNodeId(gMaxDualNodeId - gMinDualNodeId + 1, -GO_ONE);
301 Array<GlobalOrdinal> local_dualNodeId2primalAggId(gMaxDualNodeId - gMinDualNodeId + 1, -GO_ONE);
303 Array<GlobalOrdinal> dualDofId2primalDofId(primalInterfaceDofRowMap->getGlobalNumElements(), -GO_ONE);
304 Array<GlobalOrdinal> local_dualDofId2primalDofId(primalInterfaceDofRowMap->getGlobalNumElements(), -GO_ONE);
307 const size_t numMyPrimalInterfaceDOFs = primalInterfaceDofRowMap->getLocalNumElements();
308 for (
size_t r = 0; r < numMyPrimalInterfaceDOFs; r += numDofsPerPrimalNode) {
309 GlobalOrdinal gPrimalRowId = primalInterfaceDofRowMap->getGlobalElement(r);
311 if (A01->getRowMap()->isNodeGlobalElement(gPrimalRowId))
313 const LocalOrdinal lPrimalRowId = A01->getRowMap()->getLocalElement(gPrimalRowId);
315 const LocalOrdinal lPrimalNodeId = lPrimalRowId / numDofsPerPrimalNode;
316 const LocalOrdinal primalAggId = primalVertex2AggId[lPrimalNodeId];
317 const GlobalOrdinal gDualDofId = A01->getDomainMap()->getGlobalElement(r);
322 "PROC: " << myRank <<
" gDualNodeId " << gDualNodeId
323 <<
" is already connected to primal nodeId "
324 << local_dualNodeId2primalNodeId[gDualNodeId - gMinDualNodeId]
325 <<
". This shouldn't be. A possible reason might be: "
326 "Check if parallel distribution of primalInterfaceDofRowMap corresponds "
327 "to the parallel distribution of subblock matrix A01.");
329 local_dualNodeId2primalNodeId[gDualNodeId - gMinDualNodeId] = gPrimalNodeId;
330 local_dualNodeId2primalAggId[gDualNodeId - gMinDualNodeId] = primalAggId;
333 const int dualNodeId2primalNodeIdSize = Teuchos::as<int>(local_dualNodeId2primalNodeId.
size());
335 &local_dualNodeId2primalNodeId[0], &dualNodeId2primalNodeId[0]);
337 &local_dualNodeId2primalAggId[0], &dualNodeId2primalAggId[0]);
342 std::vector<GlobalOrdinal> dualNodes;
343 for (
size_t r = 0; r < A01->getDomainMap()->getLocalNumElements(); r++) {
347 GlobalOrdinal gDualDofId = A01->getDomainMap()->getGlobalElement(r);
349 dualNodes.push_back(gDualNodeId);
353 dualNodes.erase(std::unique(dualNodes.begin(), dualNodes.end()), dualNodes.end());
361 dualAggregates->setObjectLabel(
"UC (dual variables)");
369 std::map<GlobalOrdinal, LocalOrdinal> primalAggId2localDualAggId;
370 for (
size_t lDualNodeID = 0; lDualNodeID < dualNodeMap->getLocalNumElements(); ++lDualNodeID) {
371 const GlobalOrdinal gDualNodeId = dualNodeMap->getGlobalElement(lDualNodeID);
372 const GlobalOrdinal primalAggId = dualNodeId2primalAggId[gDualNodeId - gMinDualNodeId];
373 if (primalAggId2localDualAggId.count(primalAggId) == 0)
374 primalAggId2localDualAggId[primalAggId] = nLocalAggregates++;
375 dualVertex2AggId[lDualNodeID] = primalAggId2localDualAggId[primalAggId];
376 dualProcWinner[lDualNodeID] = myRank;
386 const size_t numMyDualNodes = dualNodeMap->getLocalNumElements();
387 for (
size_t lDualNodeID = 0; lDualNodeID < numMyDualNodes; ++lDualNodeID) {
388 for (
LocalOrdinal dof = 0; dof < numDofsPerDualNode; ++dof) {
389 rowTranslation->push_back(lDualNodeID);
390 colTranslation->push_back(lDualNodeID);
397 A01->getDomainMap(), A01->getDomainMap(), A01->getDomainMap(),
398 fullblocksize, dualDofOffset, blockid, nStridedOffset, stridedblocksize));
400 dualAggregates->SetNumAggregates(nLocalAggregates);
401 dualAggregates->AggregatesCrossProcessors(primalAggregates->AggregatesCrossProcessors());
403 if (dualAggregates->AggregatesCrossProcessors())
404 GetOStream(
Runtime1) <<
"Interface aggregates cross processor boundaries." << std::endl;
406 GetOStream(
Runtime1) <<
"Interface aggregates do not cross processor boundaries." << std::endl;
408 currentLevel.Set(
"Aggregates", dualAggregates,
this);
409 currentLevel.Set(
"UnAmalgamationInfo", dualAmalgamationInfo,
this);
MueLu::DefaultLocalOrdinal LocalOrdinal
const RCP< LOVector > & GetProcWinner() const
Returns constant vector that maps local node IDs to owning processor IDs.
Container class for aggregation information.
void BuildBasedOnNodeMapping(const std::string &prefix, Level ¤tLevel) const
Build dual aggregates based on a given dual-to-primal node mapping.
T & get(const std::string &name, T def_value)
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
ParameterList & set(std::string const &name, T &&value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
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.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Class that holds all level-specific information.
void Build(Level ¤tLevel) const override
Build aggregates.
virtual void setObjectLabel(const std::string &objectLabel)
const RCP< LOMultiVector > & GetVertex2AggId() const
Returns constant vector that maps local node IDs to local aggregates IDs.
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.
KOKKOS_INLINE_FUNCTION void AggregatesCrossProcessors(const bool &flag)
Record whether aggregates include DOFs from other processes.
int GetLevelID() const
Return level number.
Exception throws to report errors in the internal logical of the program.
Print all warning messages.
#define TEUCHOS_ASSERT(assertion_test)
Description of what is happening (more verbose)
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.
std::string description() const
Return a simple one-line description of this object.
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
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
void SetNumAggregates(LO nAggregates)
Set number of local aggregates on current processor.
Exception throws to report invalid user entry.