MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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_Level.hpp"
57 #include "MueLu_Monitor.hpp"
58 
60 
61 namespace MueLu {
62 
63 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
65  RCP<ParameterList> validParamList = rcp(new ParameterList());
66 
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)");
69 
70  validParamList->set<std::string>("Dual/primal mapping strategy", "vague",
71  "Strategy to represent mapping between dual and primal quantities [node-based, dof-based]");
72 
73  validParamList->set<RCP<const FactoryBase>>("DualNodeID2PrimalNodeID", Teuchos::null,
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");
75  validParamList->set<LocalOrdinal>("number of DOFs per dual node", -Teuchos::ScalarTraits<LocalOrdinal>::one(),
76  "Number of DOFs per dual node");
77 
78  validParamList->set<RCP<const FactoryBase>>("Primal interface DOF map", Teuchos::null,
79  "Generating factory of the primal DOF row map of slave side of the coupling surface");
80 
81  return validParamList;
82 } // GetValidParameterList()
83 
84 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
86  Input(currentLevel, "A"); // matrix block of dual variables
87  Input(currentLevel, "Aggregates");
88 
89  const ParameterList &pL = GetParameterList();
90  TEUCHOS_TEST_FOR_EXCEPTION(pL.get<std::string>("Dual/primal mapping strategy") == "vague", Exceptions::InvalidArgument,
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") {
93  if (currentLevel.GetLevelID() == 0) {
94  TEUCHOS_TEST_FOR_EXCEPTION(!currentLevel.IsAvailable("DualNodeID2PrimalNodeID", NoFactory::get()),
95  Exceptions::RuntimeError, "DualNodeID2PrimalNodeID was not provided by the user on level 0!");
96 
97  currentLevel.DeclareInput("DualNodeID2PrimalNodeID", NoFactory::get(), this);
98  } else {
99  Input(currentLevel, "DualNodeID2PrimalNodeID");
100  }
101  } else if (pL.get<std::string>("Dual/primal mapping strategy") == "dof-based") {
102  if (currentLevel.GetLevelID() == 0)
103  currentLevel.DeclareInput("Primal interface DOF map", NoFactory::get(), this);
104  else
105  Input(currentLevel, "Primal interface DOF map");
106  } else
107  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::InvalidArgument, "Unknown strategy for dual/primal mapping.")
108 
109 } // DeclareInput
110 
111 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
113  const std::string prefix = "MueLu::InterfaceAggregationFactory::Build: ";
114 
115  FactoryMonitor m(*this, "Build", currentLevel);
116 
117  // Call a specialized build routine based on the format of user-given input
118  const ParameterList &pL = GetParameterList();
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);
124  else
126  "MueLu::InterfaceAggregationFactory::Builld(): Unknown strategy for dual/primal mapping. Set a valid value for the parameter \"" << parameterName << "\".")
127 }
128 
129 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
131  Level &currentLevel) const {
132  using Dual2Primal_type = std::map<LocalOrdinal, LocalOrdinal>;
133 
134  const ParameterList &pL = GetParameterList();
135 
136  RCP<const Matrix> A = Get<RCP<Matrix>>(currentLevel, "A");
137  const LocalOrdinal numDofsPerDualNode = pL.get<LocalOrdinal>("number of DOFs per dual node");
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.");
140 
141  RCP<const Aggregates> primalAggregates = Get<RCP<Aggregates>>(currentLevel, "Aggregates");
142  ArrayRCP<const LocalOrdinal> primalVertex2AggId = primalAggregates->GetVertex2AggId()->getData(0);
143 
144  // Get the user-prescribed mapping of dual to primal node IDs
145  RCP<Dual2Primal_type> mapNodesDualToPrimal;
146  if (currentLevel.GetLevelID() == 0)
147  mapNodesDualToPrimal = currentLevel.Get<RCP<Dual2Primal_type>>("DualNodeID2PrimalNodeID", NoFactory::get());
148  else
149  mapNodesDualToPrimal = Get<RCP<Dual2Primal_type>>(currentLevel, "DualNodeID2PrimalNodeID");
150 
151  RCP<const Map> operatorRangeMap = A->getRangeMap();
152  const size_t myRank = operatorRangeMap->getComm()->getRank();
153 
154  LocalOrdinal globalNumDualNodes = operatorRangeMap->getGlobalNumElements() / numDofsPerDualNode;
155  LocalOrdinal localNumDualNodes = operatorRangeMap->getLocalNumElements() / numDofsPerDualNode;
156 
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.");
159 
160  RCP<const Map> dualNodeMap = Teuchos::null;
161  if (numDofsPerDualNode == 1)
162  dualNodeMap = operatorRangeMap;
163  else {
164  GlobalOrdinal indexBase = operatorRangeMap->getIndexBase();
165  auto comm = operatorRangeMap->getComm();
166  std::vector<GlobalOrdinal> myDualNodes = {};
167 
168  for (size_t i = 0; i < operatorRangeMap->getLocalNumElements(); i += numDofsPerDualNode)
169  myDualNodes.push_back((operatorRangeMap->getGlobalElement(i) - indexBase) / numDofsPerDualNode + indexBase);
170 
171  dualNodeMap = MapFactory::Build(operatorRangeMap->lib(), globalNumDualNodes, myDualNodes, indexBase, comm);
172  }
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.");
175 
176  RCP<Aggregates> dualAggregates = rcp(new Aggregates(dualNodeMap));
177  dualAggregates->setObjectLabel("InterfaceAggregation");
178 
179  // Copy setting from primal aggregates, as we copy the interface part of primal aggregates anyways
180  dualAggregates->AggregatesCrossProcessors(primalAggregates->AggregatesCrossProcessors());
181 
182  ArrayRCP<LocalOrdinal> dualVertex2AggId = dualAggregates->GetVertex2AggId()->getDataNonConst(0);
183  ArrayRCP<LocalOrdinal> dualProcWinner = dualAggregates->GetProcWinner()->getDataNonConst(0);
184 
185  RCP<Dual2Primal_type> coarseMapNodesDualToPrimal = rcp(new Dual2Primal_type());
186  RCP<Dual2Primal_type> coarseMapNodesPrimalToDual = rcp(new Dual2Primal_type());
187 
188  LocalOrdinal numLocalDualAggregates = 0;
189 
190  /* Loop over the local dual nodes and
191  *
192  * - assign dual nodes to dual aggregates
193  * - recursively coarsen the dual-to-primal node mapping
194  */
197  for (LocalOrdinal localDualNodeID = 0; localDualNodeID < localNumDualNodes; ++localDualNodeID) {
198  // Get local ID of the primal node associated to the current dual node
199  localPrimalNodeID = (*mapNodesDualToPrimal)[localDualNodeID];
200 
201  // Find the primal aggregate that owns the current primal node
202  currentPrimalAggId = primalVertex2AggId[localPrimalNodeID];
203 
204  // Test if the current primal aggregate has no associated dual aggregate, yet.
205  // Create new dual aggregate, if necessary.
206  if (coarseMapNodesPrimalToDual->count(currentPrimalAggId) == 0) {
207  // Associate a new dual aggregate w/ the current primal aggregate
208  (*coarseMapNodesPrimalToDual)[currentPrimalAggId] = numLocalDualAggregates;
209  (*coarseMapNodesDualToPrimal)[numLocalDualAggregates] = currentPrimalAggId;
210  ++numLocalDualAggregates;
211  }
212 
213  // Fill the dual aggregate
214  dualVertex2AggId[localDualNodeID] = (*coarseMapNodesPrimalToDual)[currentPrimalAggId];
215  dualProcWinner[localDualNodeID] = myRank;
216  }
217 
218  // Store dual aggregeate data as well as coarsening information
219  dualAggregates->SetNumAggregates(numLocalDualAggregates);
220  Set(currentLevel, "Aggregates", dualAggregates);
221  Set(currentLevel, "CoarseDualNodeID2PrimalNodeID", coarseMapNodesDualToPrimal);
222  GetOStream(Statistics1) << dualAggregates->description() << std::endl;
223 } // BuildBasedOnNodeMapping
224 
225 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
227  const std::string &prefix, Level &currentLevel) const {
230 
231  // filled with striding information from A01
232  LocalOrdinal numDofsPerDualNode = 0;
233  LocalOrdinal numDofsPerPrimalNode = 0;
234 
235  // Grab the off-diagonal block (0,1) from the global blocked operator
236  RCP<const Matrix> A01 = Get<RCP<Matrix>>(currentLevel, "A");
237  RCP<const Aggregates> primalAggregates = Get<RCP<Aggregates>>(currentLevel, "Aggregates");
238  ArrayRCP<const LocalOrdinal> primalVertex2AggId = primalAggregates->GetVertex2AggId()->getData(0);
239 
240  auto comm = A01->getRowMap()->getComm();
241  const int myRank = comm->getRank();
242 
243  RCP<const Map> primalInterfaceDofRowMap = Teuchos::null;
244  if (currentLevel.GetLevelID() == 0) {
245  // Use NoFactory, since the fine level asks for user data
246  primalInterfaceDofRowMap = currentLevel.Get<RCP<const Map>>("Primal interface DOF map", NoFactory::get());
247  } else {
248  primalInterfaceDofRowMap = Get<RCP<const Map>>(currentLevel, "Primal interface DOF map");
249  }
250  TEUCHOS_ASSERT(!primalInterfaceDofRowMap.is_null());
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());
256 
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;
261  }
262  }
263 
264  TEUCHOS_TEST_FOR_EXCEPTION(numDofsPerPrimalNode == 0, Exceptions::RuntimeError,
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.");
266  TEUCHOS_TEST_FOR_EXCEPTION(numDofsPerDualNode == 0, Exceptions::RuntimeError,
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.");
268 
269  /* Determine block information for primal block
270  *
271  * primalDofOffset: global offset of primal DOF GIDs (usually is zero (default))
272  * primalBlockDim: block dim for fixed size blocks
273  * - is 2 or 3 (for 2d or 3d problems) on the finest level (# displacement dofs per node)
274  * - is 3 or 6 (for 2d or 3d problems) on coarser levels (# nullspace vectors)
275  */
276  GlobalOrdinal primalDofOffset = GO_ZERO;
277  LocalOrdinal primalBlockDim = numDofsPerPrimalNode;
278 
279  /* Determine block information for Lagrange multipliers
280  *
281  * dualDofOffset: usually > zero (set by domainOffset for Ptent11Fact)
282  * dualBlockDim:
283  * - is primalBlockDim (for 2d or 3d problems) on the finest level (1 Lagrange multiplier per
284  * displacement dof)
285  * - is 2 or 3 (for 2d or 3d problems) on coarser levels (same as on finest level, whereas there
286  * are 3 or 6 displacement dofs per node)
287  */
288  GlobalOrdinal dualDofOffset = A01->getRowMap()->getMaxAllGlobalIndex() + 1;
289  LocalOrdinal dualBlockDim = numDofsPerDualNode;
290  // Generate global replicated mapping "lagrNodeId -> dispNodeId"
291  RCP<const Map> dualDofMap = A01->getDomainMap();
293  dualDofMap->getMaxAllGlobalIndex(), dualBlockDim, dualDofOffset, dualDofMap->getIndexBase());
295  dualDofMap->getMinAllGlobalIndex(), dualBlockDim, dualDofOffset, dualDofMap->getIndexBase());
296 
297  GetOStream(Runtime1) << " Dual DOF map: index base = " << dualDofMap->getIndexBase()
298  << ", block dim = " << dualBlockDim
299  << ", gid offset = " << dualDofOffset
300  << std::endl;
301 
302  GetOStream(Runtime1) << " [primal / dual] DOFs per node = [" << numDofsPerPrimalNode
303  << "/" << numDofsPerDualNode << "]" << std::endl;
304 
305  // Generate locally replicated vector for mapping dual node IDs to primal node IDs
306  Array<GlobalOrdinal> dualNodeId2primalNodeId(gMaxDualNodeId - gMinDualNodeId + 1, -GO_ONE);
307  Array<GlobalOrdinal> local_dualNodeId2primalNodeId(gMaxDualNodeId - gMinDualNodeId + 1, -GO_ONE);
308 
309  // Generate locally replicated vector for mapping dual node IDs to primal aggregate ID
310  Array<GlobalOrdinal> dualNodeId2primalAggId(gMaxDualNodeId - gMinDualNodeId + 1, -GO_ONE);
311  Array<GlobalOrdinal> local_dualNodeId2primalAggId(gMaxDualNodeId - gMinDualNodeId + 1, -GO_ONE);
312 
313  Array<GlobalOrdinal> dualDofId2primalDofId(primalInterfaceDofRowMap->getGlobalNumElements(), -GO_ONE);
314  Array<GlobalOrdinal> local_dualDofId2primalDofId(primalInterfaceDofRowMap->getGlobalNumElements(), -GO_ONE);
315 
316  // Fill mapping of Lagrange Node IDs to displacement aggregate IDs
317  const size_t numMyPrimalInterfaceDOFs = primalInterfaceDofRowMap->getLocalNumElements();
318  for (size_t r = 0; r < numMyPrimalInterfaceDOFs; r += numDofsPerPrimalNode) {
319  GlobalOrdinal gPrimalRowId = primalInterfaceDofRowMap->getGlobalElement(r);
320 
321  if (A01->getRowMap()->isNodeGlobalElement(gPrimalRowId)) // Remove this if?
322  {
323  const LocalOrdinal lPrimalRowId = A01->getRowMap()->getLocalElement(gPrimalRowId);
324  const GlobalOrdinal gPrimalNodeId = AmalgamationFactory::DOFGid2NodeId(gPrimalRowId, primalBlockDim, primalDofOffset, primalInterfaceDofRowMap->getIndexBase());
325  const LocalOrdinal lPrimalNodeId = lPrimalRowId / numDofsPerPrimalNode;
326  const LocalOrdinal primalAggId = primalVertex2AggId[lPrimalNodeId];
327  const GlobalOrdinal gDualDofId = A01->getDomainMap()->getGlobalElement(r);
328  const GlobalOrdinal gDualNodeId = AmalgamationFactory::DOFGid2NodeId(gDualDofId, dualBlockDim, dualDofOffset, 0);
329 
330  TEUCHOS_TEST_FOR_EXCEPTION(local_dualNodeId2primalNodeId[gDualNodeId - gMinDualNodeId] != -GO_ONE,
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.");
338 
339  local_dualNodeId2primalNodeId[gDualNodeId - gMinDualNodeId] = gPrimalNodeId;
340  local_dualNodeId2primalAggId[gDualNodeId - gMinDualNodeId] = primalAggId;
341  }
342  }
343  const int dualNodeId2primalNodeIdSize = Teuchos::as<int>(local_dualNodeId2primalNodeId.size());
344  Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, dualNodeId2primalNodeIdSize,
345  &local_dualNodeId2primalNodeId[0], &dualNodeId2primalNodeId[0]);
346  Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, dualNodeId2primalNodeIdSize,
347  &local_dualNodeId2primalAggId[0], &dualNodeId2primalAggId[0]);
348 
349  // build node map for dual variables
350  // generate "artificial nodes" for lagrange multipliers
351  // the node map is also used for defining the Aggregates for the lagrange multipliers
352  std::vector<GlobalOrdinal> dualNodes;
353  for (size_t r = 0; r < A01->getDomainMap()->getLocalNumElements(); r++) {
354  // determine global Lagrange multiplier row Dof
355  // generate a node id using the grid, lagr_blockdim and lagr_offset // todo make sure, that
356  // nodeId is unique and does not interfer with the displacement nodes
357  GlobalOrdinal gDualDofId = A01->getDomainMap()->getGlobalElement(r);
358  GlobalOrdinal gDualNodeId = AmalgamationFactory::DOFGid2NodeId(gDualDofId, dualBlockDim, dualDofOffset, 0);
359  dualNodes.push_back(gDualNodeId);
360  }
361 
362  // remove all duplicates
363  dualNodes.erase(std::unique(dualNodes.begin(), dualNodes.end()), dualNodes.end());
364 
365  // define node map for Lagrange multipliers
366  Teuchos::RCP<const Map> dualNodeMap = MapFactory::Build(A01->getRowMap()->lib(),
367  Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(), dualNodes, A01->getRowMap()->getIndexBase(), comm);
368 
369  // Build aggregates using the lagrange multiplier node map
370  Teuchos::RCP<Aggregates> dualAggregates = Teuchos::rcp(new Aggregates(dualNodeMap));
371  dualAggregates->setObjectLabel("UC (dual variables)");
372 
373  // extract aggregate data structures to fill
374  Teuchos::ArrayRCP<LocalOrdinal> dualVertex2AggId = dualAggregates->GetVertex2AggId()->getDataNonConst(0);
375  Teuchos::ArrayRCP<LocalOrdinal> dualProcWinner = dualAggregates->GetProcWinner()->getDataNonConst(0);
376 
377  // loop over local lagrange multiplier node ids
378  LocalOrdinal nLocalAggregates = 0;
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;
387  }
388 
389  const LocalOrdinal fullblocksize = numDofsPerDualNode;
390  const LocalOrdinal blockid = -1;
391  const LocalOrdinal nStridedOffset = 0;
392  const LocalOrdinal stridedblocksize = fullblocksize;
393 
394  RCP<Array<LO>> rowTranslation = rcp(new Array<LO>());
395  RCP<Array<LO>> colTranslation = rcp(new Array<LO>());
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);
401  }
402  }
403 
404  TEUCHOS_ASSERT(A01->isFillComplete());
405 
406  RCP<AmalgamationInfo> dualAmalgamationInfo = rcp(new AmalgamationInfo(rowTranslation, colTranslation,
407  A01->getDomainMap(), A01->getDomainMap(), A01->getDomainMap(),
408  fullblocksize, dualDofOffset, blockid, nStridedOffset, stridedblocksize));
409 
410  dualAggregates->SetNumAggregates(nLocalAggregates);
411  dualAggregates->AggregatesCrossProcessors(primalAggregates->AggregatesCrossProcessors());
412 
413  if (dualAggregates->AggregatesCrossProcessors())
414  GetOStream(Runtime1) << "Interface aggregates cross processor boundaries." << std::endl;
415  else
416  GetOStream(Runtime1) << "Interface aggregates do not cross processor boundaries." << std::endl;
417 
418  currentLevel.Set("Aggregates", dualAggregates, this);
419  currentLevel.Set("UnAmalgamationInfo", dualAmalgamationInfo, this);
420 
421 } // BuildBasedOnPrimalInterfaceDofMap
422 
423 } // namespace MueLu
424 
425 #endif /* MUELU_INTERFACEAGGREGATIONFACTORY_DEF_HPP_ */
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 &currentLevel) 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)
Print more statistics.
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.
Definition: MueLu_Level.hpp:99
void Build(Level &currentLevel) 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 &currentLevel) const override
Specifies the data that this class needs, and the factories that generate that data.
size_type size() const
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.
Definition: MueLu_Level.cpp:76
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 &currentLevel) 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&#39;s value has been saved.
void SetNumAggregates(LO nAggregates)
Set number of local aggregates on current processor.
Exception throws to report invalid user entry.
bool is_null() const