MueLu Version of the Day
Loading...
Searching...
No Matches
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
62namespace MueLu
63{
64
65template <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
87template <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
122template <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
141template <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->getLocalNumElements() / 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->getLocalNumElements(); 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->getLocalNumElements()),
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
241template <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->getLocalNumElements();
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()->getLocalNumElements(); 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->getLocalNumElements(); ++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->getLocalNumElements();
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
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Container class for aggregation information.
static const GlobalOrdinal DOFGid2NodeId(GlobalOrdinal gid, LocalOrdinal blockSize, const GlobalOrdinal offset, const GlobalOrdinal indexBase)
Translate global (row/column) id to global amalgamation block id.
minimal container class for storing amalgamation information
Exception throws to report invalid user entry.
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
void Build(Level &currentLevel) const override
Build aggregates.
void BuildBasedOnNodeMapping(const std::string &prefix, Level &currentLevel) const
Build dual aggregates based on a given dual-to-primal node mapping.
void DeclareInput(Level &currentLevel) const override
Specifies the data that this class needs, and the factories that generate that data.
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.
RCP< const ParameterList > GetValidParameterList() const override
Input.
Class that holds all level-specific information.
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
int GetLevelID() const
Return level number.
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
static const NoFactory * get()
Namespace for MueLu classes and methods.
@ Warnings
Print all warning messages.
@ Statistics1
Print more statistics.
@ Runtime1
Description of what is happening (more verbose)