MueLu  Version of the Day
MueLu_InterfaceAggregationFactory_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // MueLu: A package for multigrid based preconditioning
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact
39 // Jonathan Hu (jhu@sandia.gov)
40 // Andrey Prokopenko (aprokop@sandia.gov)
41 // Ray Tuminaro (rstumin@sandia.gov)
42 //
43 // ***********************************************************************
44 //
45 // @HEADER
46 #ifndef MUELU_INTERFACEAGGREGATIONFACTORY_DEF_HPP_
47 #define MUELU_INTERFACEAGGREGATIONFACTORY_DEF_HPP_
48 
50 
51 #include "MueLu_Level.hpp"
52 #include "MueLu_Aggregates.hpp"
53 #include "MueLu_Monitor.hpp"
54 #include "Xpetra_MapFactory.hpp"
55 
56 namespace MueLu
57 {
58 
59 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
61 {
62  RCP<ParameterList> validParamList = rcp(new ParameterList());
63 
64  validParamList->set<RCP<const FactoryBase>>("A", Teuchos::null, "Generating factory of A");
65  validParamList->set<RCP<const FactoryBase>>("Aggregates", Teuchos::null, "Generating factory of the Aggregates (for block 0,0)");
66  validParamList->set<RCP<const FactoryBase>>("DualNodeID2PrimalNodeID", Teuchos::null, "Generating factory of the DualNodeID2PrimalNodeID map");
67 
68  validParamList->set<LocalOrdinal>("number of DOFs per dual node", Teuchos::ScalarTraits<LocalOrdinal>::one(), "Number of DOFs per dual node");
69 
70  return validParamList;
71 }
72 
73 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
75 {
76  Input(currentLevel, "A");
77  Input(currentLevel, "Aggregates");
78 
79  if (currentLevel.GetLevelID() == 0)
80  {
81  if(currentLevel.IsAvailable("DualNodeID2PrimalNodeID", NoFactory::get())) {
82  currentLevel.DeclareInput("DualNodeID2PrimalNodeID", NoFactory::get(), this);
83  } else {
84  TEUCHOS_TEST_FOR_EXCEPTION(currentLevel.IsAvailable("DualNodeID2PrimalNodeID", NoFactory::get()),
86  "DualNodeID2PrimalNodeID was not provided by the user on level 0!");
87  }
88  }
89  else
90  {
91  Input(currentLevel, "DualNodeID2PrimalNodeID");
92  }
93 }
94 
95 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
97 {
98  using Map = Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>;
99  using MapFactory = Xpetra::MapFactory<LocalOrdinal, GlobalOrdinal, Node>;
101  using Dual2Primal_type = std::map<LocalOrdinal, LocalOrdinal>;
102 
103  const char prefix[] = "MueLu::InterfaceAggregationFactory::Build: ";
104  const ParameterList &pL = GetParameterList();
105 
106  FactoryMonitor m(*this, "Build", currentLevel);
107 
108  RCP<Matrix> A = Get<RCP<Matrix>>(currentLevel, "A");
109  RCP<Aggregates> aggs00 = Get<RCP<Aggregates>>(currentLevel, "Aggregates");
110  ArrayRCP<LocalOrdinal> vertex2AggIdInput = aggs00->GetVertex2AggId()->getDataNonConst(0);
111 
112  ArrayRCP<LocalOrdinal> procWinnerInput = aggs00->GetProcWinner()->getDataNonConst(0);
113 
114  RCP<Dual2Primal_type> lagr2dof;
115 
116  if (currentLevel.GetLevelID() == 0)
117  lagr2dof = currentLevel.Get<RCP<Dual2Primal_type>>("DualNodeID2PrimalNodeID", NoFactory::get());
118  else
119  lagr2dof = Get<RCP<Dual2Primal_type>>(currentLevel, "DualNodeID2PrimalNodeID");
120 
121  const LocalOrdinal nDOFPerDualNode = pL.get<LocalOrdinal>("number of DOFs per dual node");
122 
123  RCP<const Map> aggRangeMap = A->getRangeMap();
124  const size_t myRank = aggRangeMap->getComm()->getRank();
125 
126  LocalOrdinal globalNumDualNodes = aggRangeMap->getGlobalNumElements() / nDOFPerDualNode;
127  LocalOrdinal numDualNodes = aggRangeMap->getNodeNumElements() / nDOFPerDualNode;
128 
129  TEUCHOS_TEST_FOR_EXCEPTION(numDualNodes != LocalOrdinal(lagr2dof->size()), std::runtime_error, prefix << " MueLu requires the range map and the DualNodeID2PrimalNodeID map to be compatible.");
130 
131  RCP<const Map> aggVertexMap;
132 
133  if (nDOFPerDualNode == 1)
134  aggVertexMap = aggRangeMap;
135  else
136  {
137  GlobalOrdinal indexBase = aggRangeMap->getIndexBase();
138  auto comm = aggRangeMap->getComm();
139  std::vector<GlobalOrdinal> myDualNodes = {};
140 
141  for (size_t i = 0; i < aggRangeMap->getNodeNumElements(); i += nDOFPerDualNode)
142  myDualNodes.push_back((aggRangeMap->getGlobalElement(i) - indexBase) / nDOFPerDualNode + indexBase);
143 
144  aggVertexMap = MapFactory::Build(aggRangeMap->lib(), globalNumDualNodes, myDualNodes, indexBase, comm);
145  }
146 
147  RCP<Aggregates> aggregates = rcp(new Aggregates(aggVertexMap));
148  aggregates->setObjectLabel("IA");
149  aggregates->AggregatesCrossProcessors(false);
150 
151  ArrayRCP<LocalOrdinal> vertex2AggId = aggregates->GetVertex2AggId()->getDataNonConst(0);
152  ArrayRCP<LocalOrdinal> procWinner = aggregates->GetProcWinner()->getDataNonConst(0);
153 
154  RCP<Dual2Primal_type> coarseLagr2Dof = rcp(new Dual2Primal_type());
155  RCP<Dual2Primal_type> coarseDof2Lagr = rcp(new Dual2Primal_type());
156 
157  LocalOrdinal numAggId = 0;
158 
159  // Loop over the local "nodes" of the block 11
160  for (LocalOrdinal localDualNodeID = 0; localDualNodeID < numDualNodes; ++localDualNodeID)
161  {
162  // Get local ID of the primal node associated to the current dual node
163  LocalOrdinal localPrimalNodeID = (*lagr2dof)[localDualNodeID];
164 
165  // Find the aggregate of block 00 associated with the current primal node
166  LocalOrdinal currentAggIdInput = vertex2AggIdInput[localPrimalNodeID];
167 
168  // Test if the current aggregate of block 00 has no associated aggregate of block 11
169  if (coarseDof2Lagr->count(currentAggIdInput) == 0)
170  {
171  // Associate a new aggregate of block 11 to the current aggregate of block 00
172  (*coarseDof2Lagr)[currentAggIdInput] = numAggId;
173  (*coarseLagr2Dof)[numAggId] = currentAggIdInput;
174  ++numAggId;
175  }
176 
177  // Fill the block 11 aggregate information
178  vertex2AggId[localDualNodeID] = (*coarseDof2Lagr)[currentAggIdInput];
179  procWinner[localDualNodeID] = myRank;
180  }
181 
182  aggregates->SetNumAggregates(numAggId);
183  Set(currentLevel, "Aggregates", aggregates);
184  Set(currentLevel, "CoarseDualNodeID2PrimalNodeID", coarseLagr2Dof);
185  GetOStream(Statistics1) << aggregates->description() << std::endl;
186 }
187 } //namespace MueLu
188 
189 #endif /* MUELU_INTERFACEAGGREGATIONFACTORY_DEF_HPP_ */
MueLu::DefaultLocalOrdinal LocalOrdinal
Container class for aggregation information.
Timer to be used in factories. Similar to Monitor but with additional timers.
Print more statistics.
Namespace for MueLu classes and methods.
static const NoFactory * get()
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:76
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
void Build(Level &currentLevel) const override
Build aggregates.
void DeclareInput(Level &currentLevel) const override
Input.
RCP< const ParameterList > GetValidParameterList() const override
Return a const parameter list of valid parameters that setParameterList() will accept.
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need&#39;s value has been saved.
Exception throws to report errors in the internal logical of the program.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()