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 
49 #include <Xpetra_Map.hpp>
50 #include <Xpetra_MapFactory.hpp>
51 #include <Xpetra_StridedMap.hpp>
52 
53 #include "MueLu_Aggregates.hpp"
54 #include "MueLu_AmalgamationFactory.hpp"
55 #include "MueLu_AmalgamationInfo.hpp"
56 #include "MueLu_FactoryManager.hpp"
57 #include "MueLu_Level.hpp"
58 #include "MueLu_Monitor.hpp"
59 
61 
62 namespace MueLu
63 {
64 
65 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
67 {
68  RCP<ParameterList> validParamList = rcp(new ParameterList());
69 
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)");
72 
73  validParamList->set<std::string>("Dual/primal mapping strategy", "vague",
74  "Strategy to represent mapping between dual and primal quantities [node-based, dof-based]");
75 
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");
80 
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");
83 
84  return validParamList;
85 } // GetValidParameterList()
86 
87 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
89 {
90  Input(currentLevel, "A"); // matrix block of dual variables
91  Input(currentLevel, "Aggregates");
92 
93  const ParameterList &pL = GetParameterList();
94  TEUCHOS_TEST_FOR_EXCEPTION(pL.get<std::string>("Dual/primal mapping strategy")=="vague", Exceptions::InvalidArgument,
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")
97  {
98  if (currentLevel.GetLevelID() == 0)
99  {
100  TEUCHOS_TEST_FOR_EXCEPTION(!currentLevel.IsAvailable("DualNodeID2PrimalNodeID", NoFactory::get()),
101  Exceptions::RuntimeError, "DualNodeID2PrimalNodeID was not provided by the user on level 0!");
102 
103  currentLevel.DeclareInput("DualNodeID2PrimalNodeID", NoFactory::get(), this);
104  }
105  else
106  {
107  Input(currentLevel, "DualNodeID2PrimalNodeID");
108  }
109  }
110  else if (pL.get<std::string>("Dual/primal mapping strategy") == "dof-based")
111  {
112  if (currentLevel.GetLevelID() == 0)
113  currentLevel.DeclareInput("Primal interface DOF map", NoFactory::get(), this);
114  else
115  Input(currentLevel, "Primal interface DOF map");
116  }
117  else
118  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::InvalidArgument, "Unknown strategy for dual/primal mapping.")
119 
120 } // DeclareInput
121 
122 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
124 {
125  const std::string prefix = "MueLu::InterfaceAggregationFactory::Build: ";
126 
127  FactoryMonitor m(*this, "Build", currentLevel);
128 
129  // Call a specialized build routine based on the format of user-given input
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);
136  else
137  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::InvalidArgument,
138  "MueLu::InterfaceAggregationFactory::Builld(): Unknown strategy for dual/primal mapping. Set a valid value for the parameter \"" << parameterName << "\".")
139 }
140 
141 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
143  Level &currentLevel) const
144 {
145  using Dual2Primal_type = std::map<LocalOrdinal, LocalOrdinal>;
146 
147  const ParameterList &pL = GetParameterList();
148 
149  RCP<const Matrix> A = Get<RCP<Matrix>>(currentLevel, "A");
150  const LocalOrdinal numDofsPerDualNode = pL.get<LocalOrdinal>("number of DOFs per dual node");
151  TEUCHOS_TEST_FOR_EXCEPTION(numDofsPerDualNode<Teuchos::ScalarTraits<LocalOrdinal>::one(), Exceptions::InvalidArgument,
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.");
153 
154  RCP<const Aggregates> primalAggregates = Get<RCP<Aggregates>>(currentLevel, "Aggregates");
155  ArrayRCP<const LocalOrdinal> primalVertex2AggId = primalAggregates->GetVertex2AggId()->getData(0);
156 
157  // Get the user-prescribed mapping of dual to primal node IDs
158  RCP<Dual2Primal_type> mapNodesDualToPrimal;
159  if (currentLevel.GetLevelID() == 0)
160  mapNodesDualToPrimal = currentLevel.Get<RCP<Dual2Primal_type>>("DualNodeID2PrimalNodeID", NoFactory::get());
161  else
162  mapNodesDualToPrimal = Get<RCP<Dual2Primal_type>>(currentLevel, "DualNodeID2PrimalNodeID");
163 
164  RCP<const Map> operatorRangeMap = A->getRangeMap();
165  const size_t myRank = operatorRangeMap->getComm()->getRank();
166 
167  LocalOrdinal globalNumDualNodes = operatorRangeMap->getGlobalNumElements() / numDofsPerDualNode;
168  LocalOrdinal localNumDualNodes = operatorRangeMap->getNodeNumElements() / numDofsPerDualNode;
169 
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.");
172 
173  RCP<const Map> dualNodeMap = Teuchos::null;
174  if (numDofsPerDualNode == 1)
175  dualNodeMap = operatorRangeMap;
176  else
177  {
178  GlobalOrdinal indexBase = operatorRangeMap->getIndexBase();
179  auto comm = operatorRangeMap->getComm();
180  std::vector<GlobalOrdinal> myDualNodes = {};
181 
182  for (size_t i = 0; i < operatorRangeMap->getNodeNumElements(); i += numDofsPerDualNode)
183  myDualNodes.push_back((operatorRangeMap->getGlobalElement(i) - indexBase) / numDofsPerDualNode + indexBase);
184 
185  dualNodeMap = MapFactory::Build(operatorRangeMap->lib(), globalNumDualNodes, myDualNodes, indexBase, comm);
186  }
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.");
189 
190  RCP<Aggregates> dualAggregates = rcp(new Aggregates(dualNodeMap));
191  dualAggregates->setObjectLabel("InterfaceAggregation");
192 
193  // Copy setting from primal aggregates, as we copy the interface part of primal aggregates anyways
194  dualAggregates->AggregatesCrossProcessors(primalAggregates->AggregatesCrossProcessors());
195 
196  ArrayRCP<LocalOrdinal> dualVertex2AggId = dualAggregates->GetVertex2AggId()->getDataNonConst(0);
197  ArrayRCP<LocalOrdinal> dualProcWinner = dualAggregates->GetProcWinner()->getDataNonConst(0);
198 
199  RCP<Dual2Primal_type> coarseMapNodesDualToPrimal = rcp(new Dual2Primal_type());
200  RCP<Dual2Primal_type> coarseMapNodesPrimalToDual = rcp(new Dual2Primal_type());
201 
202  LocalOrdinal numLocalDualAggregates = 0;
203 
204  /* Loop over the local dual nodes and
205  *
206  * - assign dual nodes to dual aggregates
207  * - recursively coarsen the dual-to-primal node mapping
208  */
209  LocalOrdinal localPrimalNodeID = - Teuchos::ScalarTraits<LocalOrdinal>::one();
210  LocalOrdinal currentPrimalAggId = - Teuchos::ScalarTraits<LocalOrdinal>::one();
211  for (LocalOrdinal localDualNodeID = 0; localDualNodeID < localNumDualNodes; ++localDualNodeID)
212  {
213  // Get local ID of the primal node associated to the current dual node
214  localPrimalNodeID = (*mapNodesDualToPrimal)[localDualNodeID];
215 
216  // Find the primal aggregate that owns the current primal node
217  currentPrimalAggId = primalVertex2AggId[localPrimalNodeID];
218 
219  // Test if the current primal aggregate has no associated dual aggregate, yet.
220  // Create new dual aggregate, if necessary.
221  if (coarseMapNodesPrimalToDual->count(currentPrimalAggId) == 0)
222  {
223  // Associate a new dual aggregate w/ the current primal aggregate
224  (*coarseMapNodesPrimalToDual)[currentPrimalAggId] = numLocalDualAggregates;
225  (*coarseMapNodesDualToPrimal)[numLocalDualAggregates] = currentPrimalAggId;
226  ++numLocalDualAggregates;
227  }
228 
229  // Fill the dual aggregate
230  dualVertex2AggId[localDualNodeID] = (*coarseMapNodesPrimalToDual)[currentPrimalAggId];
231  dualProcWinner[localDualNodeID] = myRank;
232  }
233 
234  // Store dual aggregeate data as well as coarsening information
235  dualAggregates->SetNumAggregates(numLocalDualAggregates);
236  Set(currentLevel, "Aggregates", dualAggregates);
237  Set(currentLevel, "CoarseDualNodeID2PrimalNodeID", coarseMapNodesDualToPrimal);
238  GetOStream(Statistics1) << dualAggregates->description() << std::endl;
239 } // BuildBasedOnNodeMapping
240 
241 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
243  const std::string& prefix, Level &currentLevel) const
244 {
245  const GlobalOrdinal GO_ZERO = Teuchos::ScalarTraits<LocalOrdinal>::zero();
246  const GlobalOrdinal GO_ONE = Teuchos::ScalarTraits<GlobalOrdinal>::one();
247 
248  // filled with striding information from A01
249  LocalOrdinal numDofsPerDualNode = 0;
250  LocalOrdinal numDofsPerPrimalNode = 0;
251 
252  // Grab the off-diagonal block (0,1) from the global blocked operator
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);
256 
257  auto comm = A01->getRowMap()->getComm();
258  const int myRank = comm->getRank();
259 
260  RCP<const Map> primalInterfaceDofRowMap = Teuchos::null;
261  if (currentLevel.GetLevelID() == 0) {
262  // Use NoFactory, since the fine level asks for user data
263  primalInterfaceDofRowMap = currentLevel.Get<RCP<const Map>>("Primal interface DOF map", NoFactory::get());
264  } else {
265  primalInterfaceDofRowMap = Get<RCP<const Map>>(currentLevel, "Primal interface DOF map");
266  }
267  TEUCHOS_ASSERT(!primalInterfaceDofRowMap.is_null());
268 
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());
274 
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;
279  }
280  }
281 
282  TEUCHOS_TEST_FOR_EXCEPTION(numDofsPerPrimalNode==0, Exceptions::RuntimeError,
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.");
284  TEUCHOS_TEST_FOR_EXCEPTION(numDofsPerDualNode==0, Exceptions::RuntimeError,
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.");
286 
287  /* Determine block information for primal block
288  *
289  * primalDofOffset: global offset of primal DOF GIDs (usually is zero (default))
290  * primalBlockDim: block dim for fixed size blocks
291  * - is 2 or 3 (for 2d or 3d problems) on the finest level (# displacement dofs per node)
292  * - is 3 or 6 (for 2d or 3d problems) on coarser levels (# nullspace vectors)
293  */
294  GlobalOrdinal primalDofOffset = GO_ZERO;
295  LocalOrdinal primalBlockDim = numDofsPerPrimalNode;
296 
297  /* Determine block information for Lagrange multipliers
298  *
299  * dualDofOffset: usually > zero (set by domainOffset for Ptent11Fact)
300  * dualBlockDim:
301  * - is primalBlockDim (for 2d or 3d problems) on the finest level (1 Lagrange multiplier per
302  * displacement dof)
303  * - is 2 or 3 (for 2d or 3d problems) on coarser levels (same as on finest level, whereas there
304  * are 3 or 6 displacement dofs per node)
305  */
306  GlobalOrdinal dualDofOffset = A01->getColMap()->getMinAllGlobalIndex();
307  LocalOrdinal dualBlockDim = numDofsPerDualNode;
308 
309  // Generate global replicated mapping "lagrNodeId -> dispNodeId"
310  RCP<const Map> dualDofMap = A01->getDomainMap();
312  dualDofMap->getMaxAllGlobalIndex(), dualBlockDim, dualDofOffset, dualDofMap->getIndexBase());
314  dualDofMap->getMinAllGlobalIndex(), dualBlockDim, dualDofOffset, dualDofMap->getIndexBase());
315 
316  GetOStream(Runtime1) << " Dual DOF map: index base = " << dualDofMap->getIndexBase()
317  << ", block dim = " << dualBlockDim
318  << ", gid offset = " << dualDofOffset
319  << std::endl;
320 
321  GetOStream(Runtime1) << " [primal / dual] DOFs per node = [" << numDofsPerPrimalNode
322  << "/" << numDofsPerDualNode << "]" << std::endl;
323 
324  // Generate locally replicated vector for mapping dual node IDs to primal node IDs
325  Array<GlobalOrdinal> dualNodeId2primalNodeId(gMaxDualNodeId - gMinDualNodeId + 1, -GO_ONE);
326  Array<GlobalOrdinal> local_dualNodeId2primalNodeId(gMaxDualNodeId - gMinDualNodeId + 1, -GO_ONE);
327 
328  // Generate locally replicated vector for mapping dual node IDs to primal aggregate ID
329  Array<GlobalOrdinal> dualNodeId2primalAggId(gMaxDualNodeId - gMinDualNodeId + 1, -GO_ONE);
330  Array<GlobalOrdinal> local_dualNodeId2primalAggId(gMaxDualNodeId - gMinDualNodeId + 1, -GO_ONE);
331 
332  Array<GlobalOrdinal> dualDofId2primalDofId(primalInterfaceDofRowMap->getGlobalNumElements(), -GO_ONE);
333  Array<GlobalOrdinal> local_dualDofId2primalDofId(primalInterfaceDofRowMap->getGlobalNumElements(), -GO_ONE);
334 
335  // Fill mapping of Lagrange Node IDs to displacement aggregate IDs
336  const size_t numMyPrimalInterfaceDOFs = primalInterfaceDofRowMap->getNodeNumElements();
337  for (size_t r = 0; r < numMyPrimalInterfaceDOFs; r += numDofsPerPrimalNode)
338  {
339  GlobalOrdinal gPrimalRowId = primalInterfaceDofRowMap->getGlobalElement(r);
340 
341  if (A01->getRowMap()->isNodeGlobalElement(gPrimalRowId)) // Remove this if?
342  {
343  const LocalOrdinal lPrimalRowId = A01->getRowMap()->getLocalElement(gPrimalRowId);
344  const GlobalOrdinal gPrimalNodeId = AmalgamationFactory::DOFGid2NodeId(gPrimalRowId, primalBlockDim, primalDofOffset, primalInterfaceDofRowMap->getIndexBase());
345  const LocalOrdinal lPrimalNodeId = lPrimalRowId / numDofsPerPrimalNode;
346  const LocalOrdinal primalAggId = primalVertex2AggId[lPrimalNodeId];
347 
348  const GlobalOrdinal gDualDofId = A01->getColMap()->getGlobalElement(r);
349 
350  const GlobalOrdinal gDualNodeId = AmalgamationFactory::DOFGid2NodeId(gDualDofId, dualBlockDim, dualDofOffset, 0);
351 
352  if (local_dualNodeId2primalNodeId[gDualNodeId - gMinDualNodeId] == -GO_ONE) {
353  local_dualNodeId2primalNodeId[gDualNodeId - gMinDualNodeId] = gPrimalNodeId;
354  local_dualNodeId2primalAggId[gDualNodeId - gMinDualNodeId] = primalAggId;
355  } else {
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;
359  }
360 
361  }
362  }
363 
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]);
369 
370  // build node map for dual variables
371  // generate "artificial nodes" for lagrange multipliers
372  // the node map is also used for defining the Aggregates for the lagrange multipliers
373  std::vector<GlobalOrdinal> dualNodes;
374  for (size_t r = 0; r < A01->getDomainMap()->getNodeNumElements(); r++)
375  {
376  // determine global Lagrange multiplier row Dof
377  // generate a node id using the grid, lagr_blockdim and lagr_offset // todo make sure, that
378  // nodeId is unique and does not interfer with the displacement nodes
379  GlobalOrdinal gDualDofId = A01->getDomainMap()->getGlobalElement(r);
380  GlobalOrdinal gDualNodeId = AmalgamationFactory::DOFGid2NodeId(gDualDofId, dualBlockDim, dualDofOffset, 0);
381  dualNodes.push_back(gDualNodeId);
382  }
383 
384  // remove all duplicates
385  dualNodes.erase(std::unique(dualNodes.begin(), dualNodes.end()), dualNodes.end());
386 
387  // define node map for Lagrange multipliers
388  Teuchos::RCP<const Map> dualNodeMap = MapFactory::Build(A01->getRowMap()->lib(),
389  Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(), dualNodes, A01->getRowMap()->getIndexBase(), comm);
390 
391  // Build aggregates using the lagrange multiplier node map
392  Teuchos::RCP<Aggregates> dualAggregates = Teuchos::rcp(new Aggregates(dualNodeMap));
393  dualAggregates->setObjectLabel("UC (dual variables)");
394 
395  // extract aggregate data structures to fill
396  Teuchos::ArrayRCP<LocalOrdinal> dualVertex2AggId = dualAggregates->GetVertex2AggId()->getDataNonConst(0);
397  Teuchos::ArrayRCP<LocalOrdinal> dualProcWinner = dualAggregates->GetProcWinner()->getDataNonConst(0);
398 
399  // loop over local lagrange multiplier node ids
400  LocalOrdinal nLocalAggregates = 0;
401  std::map<GlobalOrdinal, LocalOrdinal> primalAggId2localDualAggId;
402  for (size_t lDualNodeID = 0; lDualNodeID < dualNodeMap->getNodeNumElements(); ++lDualNodeID)
403  {
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;
410  }
411 
412  const LocalOrdinal fullblocksize = numDofsPerDualNode;
413  const GlobalOrdinal offset = A01->getColMap()->getMinAllGlobalIndex();
414  const LocalOrdinal blockid = -1;
415  const LocalOrdinal nStridedOffset = 0;
416  const LocalOrdinal stridedblocksize = fullblocksize;
417 
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);
425  }
426  }
427 
428  TEUCHOS_ASSERT(A01->isFillComplete());
429 
430  RCP<AmalgamationInfo> dualAmalgamationInfo = rcp(new AmalgamationInfo(rowTranslation, colTranslation,
431  A01->getDomainMap(), A01->getDomainMap(), A01->getDomainMap(),
432  fullblocksize, offset, blockid, nStridedOffset, stridedblocksize));
433 
434  dualAggregates->SetNumAggregates(nLocalAggregates);
435  dualAggregates->AggregatesCrossProcessors(primalAggregates->AggregatesCrossProcessors());
436 
437  if (dualAggregates->AggregatesCrossProcessors())
438  GetOStream(Runtime1) << "Interface aggregates cross processor boundaries." << std::endl;
439  else
440  GetOStream(Runtime1) << "Interface aggregates do not cross processor boundaries." << std::endl;
441 
442  currentLevel.Set("Aggregates", dualAggregates, this);
443  currentLevel.Set("UnAmalgamationInfo", dualAmalgamationInfo, this);
444 
445 } // BuildBasedOnPrimalInterfaceDofMap
446 
447 } // namespace MueLu
448 
449 #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()
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.
Definition: MueLu_Level.cpp:76
void BuildBasedOnNodeMapping(const std::string &prefix, Level &currentLevel) const
Build dual aggregates based on a given dual-to-primal node mapping.
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
void Build(Level &currentLevel) const override
Build aggregates.
void BuildBasedOnPrimalInterfaceDofMap(const std::string &prefix, Level &currentLevel) const
Build dual aggregates based on a given interface row map of the primal and dual problem.
void DeclareInput(Level &currentLevel) 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&#39;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.