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"
63 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
67 validParamList->
set<
RCP<const FactoryBase>>(
"A", Teuchos::null,
"Generating factory of A (matrix block related to dual DOFs)");
68 validParamList->
set<
RCP<const FactoryBase>>(
"Aggregates", Teuchos::null,
"Generating factory of the Aggregates (for block 0,0)");
70 validParamList->
set<std::string>(
"Dual/primal mapping strategy",
"vague",
71 "Strategy to represent mapping between dual and primal quantities [node-based, dof-based]");
74 "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");
76 "Number of DOFs per dual node");
79 "Generating factory of the primal DOF row map of slave side of the coupling surface");
81 return validParamList;
84 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
86 Input(currentLevel,
"A");
87 Input(currentLevel,
"Aggregates");
91 "Strategy for dual/primal mapping not selected. Please select one of the available strategies.")
92 if (pL.
get<std::string>(
"Dual/primal mapping strategy") ==
"node-based") {
99 Input(currentLevel,
"DualNodeID2PrimalNodeID");
101 }
else if (pL.
get<std::string>(
"Dual/primal mapping strategy") ==
"dof-based") {
105 Input(currentLevel,
"Primal interface DOF map");
111 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
113 const std::string prefix =
"MueLu::InterfaceAggregationFactory::Build: ";
119 const std::string parameterName =
"Dual/primal mapping strategy";
120 if (pL.
get<std::string>(parameterName) ==
"node-based")
121 BuildBasedOnNodeMapping(prefix, currentLevel);
122 else if (pL.
get<std::string>(parameterName) ==
"dof-based")
123 BuildBasedOnPrimalInterfaceDofMap(prefix, currentLevel);
126 "MueLu::InterfaceAggregationFactory::Builld(): Unknown strategy for dual/primal mapping. Set a valid value for the parameter \"" << parameterName <<
"\".")
129 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
131 Level ¤tLevel)
const {
132 using Dual2Primal_type = std::map<LocalOrdinal, LocalOrdinal>;
139 "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.");
146 if (currentLevel.GetLevelID() == 0)
149 mapNodesDualToPrimal = Get<RCP<Dual2Primal_type>>(currentLevel,
"DualNodeID2PrimalNodeID");
152 const size_t myRank = operatorRangeMap->getComm()->getRank();
154 LocalOrdinal globalNumDualNodes = operatorRangeMap->getGlobalNumElements() / numDofsPerDualNode;
155 LocalOrdinal localNumDualNodes = operatorRangeMap->getLocalNumElements() / numDofsPerDualNode;
157 TEUCHOS_TEST_FOR_EXCEPTION(localNumDualNodes != Teuchos::as<LocalOrdinal>(mapNodesDualToPrimal->size()),
158 std::runtime_error, prefix <<
" MueLu requires the range map and the DualNodeID2PrimalNodeID map to be compatible.");
161 if (numDofsPerDualNode == 1)
162 dualNodeMap = operatorRangeMap;
165 auto comm = operatorRangeMap->getComm();
166 std::vector<GlobalOrdinal> myDualNodes = {};
168 for (
size_t i = 0; i < operatorRangeMap->getLocalNumElements(); i += numDofsPerDualNode)
169 myDualNodes.push_back((operatorRangeMap->getGlobalElement(i) - indexBase) / numDofsPerDualNode + indexBase);
171 dualNodeMap = MapFactory::Build(operatorRangeMap->lib(), globalNumDualNodes, myDualNodes, indexBase, comm);
173 TEUCHOS_TEST_FOR_EXCEPTION(localNumDualNodes != Teuchos::as<LocalOrdinal>(dualNodeMap->getLocalNumElements()),
174 std::runtime_error, prefix <<
" Local number of dual nodes given by user is incompatible to the dual node map.");
197 for (
LocalOrdinal localDualNodeID = 0; localDualNodeID < localNumDualNodes; ++localDualNodeID) {
199 localPrimalNodeID = (*mapNodesDualToPrimal)[localDualNodeID];
202 currentPrimalAggId = primalVertex2AggId[localPrimalNodeID];
206 if (coarseMapNodesPrimalToDual->count(currentPrimalAggId) == 0) {
208 (*coarseMapNodesPrimalToDual)[currentPrimalAggId] = numLocalDualAggregates;
209 (*coarseMapNodesDualToPrimal)[numLocalDualAggregates] = currentPrimalAggId;
210 ++numLocalDualAggregates;
214 dualVertex2AggId[localDualNodeID] = (*coarseMapNodesPrimalToDual)[currentPrimalAggId];
215 dualProcWinner[localDualNodeID] = myRank;
220 Set(currentLevel,
"Aggregates", dualAggregates);
221 Set(currentLevel,
"CoarseDualNodeID2PrimalNodeID", coarseMapNodesDualToPrimal);
225 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
227 const std::string &prefix,
Level ¤tLevel)
const {
240 auto comm = A01->getRowMap()->getComm();
241 const int myRank = comm->getRank();
244 if (currentLevel.GetLevelID() == 0) {
248 primalInterfaceDofRowMap = Get<RCP<const Map>>(currentLevel,
"Primal interface DOF map");
251 if (A01->IsView(
"stridedMaps") && rcp_dynamic_cast<
const StridedMap>(A01->getRowMap(
"stridedMaps")) != Teuchos::null) {
252 auto stridedRowMap = rcp_dynamic_cast<
const StridedMap>(A01->getRowMap(
"stridedMaps"));
253 auto stridedColMap = rcp_dynamic_cast<
const StridedMap>(A01->getColMap(
"stridedMaps"));
254 numDofsPerPrimalNode = Teuchos::as<LocalOrdinal>(stridedRowMap->getFixedBlockSize());
255 numDofsPerDualNode = Teuchos::as<LocalOrdinal>(stridedColMap->getFixedBlockSize());
257 if (numDofsPerPrimalNode != numDofsPerDualNode) {
258 GetOStream(
Warnings) <<
"InterfaceAggregation attempts to work with "
259 << numDofsPerPrimalNode <<
" primal DOFs per node and " << numDofsPerDualNode <<
" dual DOFs per node."
260 <<
"Be careful! Algorithm is not well-tested, if number of primal and dual DOFs per node differ." << std::endl;
265 "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.");
267 "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.");
288 GlobalOrdinal dualDofOffset = A01->getRowMap()->getMaxAllGlobalIndex() + 1;
293 dualDofMap->getMaxAllGlobalIndex(), dualBlockDim, dualDofOffset, dualDofMap->getIndexBase());
295 dualDofMap->getMinAllGlobalIndex(), dualBlockDim, dualDofOffset, dualDofMap->getIndexBase());
297 GetOStream(
Runtime1) <<
" Dual DOF map: index base = " << dualDofMap->getIndexBase()
298 <<
", block dim = " << dualBlockDim
299 <<
", gid offset = " << dualDofOffset
302 GetOStream(
Runtime1) <<
" [primal / dual] DOFs per node = [" << numDofsPerPrimalNode
303 <<
"/" << numDofsPerDualNode <<
"]" << std::endl;
307 Array<GlobalOrdinal> local_dualNodeId2primalNodeId(gMaxDualNodeId - gMinDualNodeId + 1, -GO_ONE);
311 Array<GlobalOrdinal> local_dualNodeId2primalAggId(gMaxDualNodeId - gMinDualNodeId + 1, -GO_ONE);
313 Array<GlobalOrdinal> dualDofId2primalDofId(primalInterfaceDofRowMap->getGlobalNumElements(), -GO_ONE);
314 Array<GlobalOrdinal> local_dualDofId2primalDofId(primalInterfaceDofRowMap->getGlobalNumElements(), -GO_ONE);
317 const size_t numMyPrimalInterfaceDOFs = primalInterfaceDofRowMap->getLocalNumElements();
318 for (
size_t r = 0; r < numMyPrimalInterfaceDOFs; r += numDofsPerPrimalNode) {
319 GlobalOrdinal gPrimalRowId = primalInterfaceDofRowMap->getGlobalElement(r);
321 if (A01->getRowMap()->isNodeGlobalElement(gPrimalRowId))
323 const LocalOrdinal lPrimalRowId = A01->getRowMap()->getLocalElement(gPrimalRowId);
325 const LocalOrdinal lPrimalNodeId = lPrimalRowId / numDofsPerPrimalNode;
326 const LocalOrdinal primalAggId = primalVertex2AggId[lPrimalNodeId];
327 const GlobalOrdinal gDualDofId = A01->getDomainMap()->getGlobalElement(r);
332 "PROC: " << myRank <<
" gDualNodeId " << gDualNodeId
333 <<
" is already connected to primal nodeId "
334 << local_dualNodeId2primalNodeId[gDualNodeId - gMinDualNodeId]
335 <<
". This shouldn't be. A possible reason might be: "
336 "Check if parallel distribution of primalInterfaceDofRowMap corresponds "
337 "to the parallel distribution of subblock matrix A01.");
339 local_dualNodeId2primalNodeId[gDualNodeId - gMinDualNodeId] = gPrimalNodeId;
340 local_dualNodeId2primalAggId[gDualNodeId - gMinDualNodeId] = primalAggId;
343 const int dualNodeId2primalNodeIdSize = Teuchos::as<int>(local_dualNodeId2primalNodeId.
size());
345 &local_dualNodeId2primalNodeId[0], &dualNodeId2primalNodeId[0]);
347 &local_dualNodeId2primalAggId[0], &dualNodeId2primalAggId[0]);
352 std::vector<GlobalOrdinal> dualNodes;
353 for (
size_t r = 0; r < A01->getDomainMap()->getLocalNumElements(); r++) {
357 GlobalOrdinal gDualDofId = A01->getDomainMap()->getGlobalElement(r);
359 dualNodes.push_back(gDualNodeId);
363 dualNodes.erase(std::unique(dualNodes.begin(), dualNodes.end()), dualNodes.end());
371 dualAggregates->setObjectLabel(
"UC (dual variables)");
379 std::map<GlobalOrdinal, LocalOrdinal> primalAggId2localDualAggId;
380 for (
size_t lDualNodeID = 0; lDualNodeID < dualNodeMap->getLocalNumElements(); ++lDualNodeID) {
381 const GlobalOrdinal gDualNodeId = dualNodeMap->getGlobalElement(lDualNodeID);
382 const GlobalOrdinal primalAggId = dualNodeId2primalAggId[gDualNodeId - gMinDualNodeId];
383 if (primalAggId2localDualAggId.count(primalAggId) == 0)
384 primalAggId2localDualAggId[primalAggId] = nLocalAggregates++;
385 dualVertex2AggId[lDualNodeID] = primalAggId2localDualAggId[primalAggId];
386 dualProcWinner[lDualNodeID] = myRank;
396 const size_t numMyDualNodes = dualNodeMap->getLocalNumElements();
397 for (
size_t lDualNodeID = 0; lDualNodeID < numMyDualNodes; ++lDualNodeID) {
398 for (
LocalOrdinal dof = 0; dof < numDofsPerDualNode; ++dof) {
399 rowTranslation->push_back(lDualNodeID);
400 colTranslation->push_back(lDualNodeID);
407 A01->getDomainMap(), A01->getDomainMap(), A01->getDomainMap(),
408 fullblocksize, dualDofOffset, blockid, nStridedOffset, stridedblocksize));
410 dualAggregates->SetNumAggregates(nLocalAggregates);
411 dualAggregates->AggregatesCrossProcessors(primalAggregates->AggregatesCrossProcessors());
413 if (dualAggregates->AggregatesCrossProcessors())
414 GetOStream(
Runtime1) <<
"Interface aggregates cross processor boundaries." << std::endl;
416 GetOStream(
Runtime1) <<
"Interface aggregates do not cross processor boundaries." << std::endl;
418 currentLevel.Set(
"Aggregates", dualAggregates,
this);
419 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)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
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.