46 #ifndef MUELU_BRICKAGGREGATIONFACTORY_DEF_HPP_
47 #define MUELU_BRICKAGGREGATIONFACTORY_DEF_HPP_
52 #include <Teuchos_CommHelpers.hpp>
58 #include <Xpetra_Map.hpp>
59 #include <Xpetra_MapFactory.hpp>
62 #include <Xpetra_MultiVectorFactory.hpp>
64 #include "MueLu_Aggregates.hpp"
68 #include "MueLu_Utilities.hpp"
69 #include "MueLu_LWGraph.hpp"
71 #include "MueLu_LWGraph.hpp"
75 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
79 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
86 #undef SET_VALID_ENTRY
90 return validParamList;
93 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
95 Input(currentLevel,
"A");
96 Input(currentLevel,
"Coordinates");
143 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
151 RCP<Matrix> A = Get<RCP<Matrix> >(currentLevel,
"A");
157 int numProcs = comm->getSize();
158 int myRank = comm->getRank();
160 int numPoints = colMap->getLocalNumElements();
162 bx_ = pL.
get<
int>(
"aggregation: brick x size");
163 by_ = pL.
get<
int>(
"aggregation: brick y size");
164 bz_ = pL.
get<
int>(
"aggregation: brick z size");
166 dirichletX_ = pL.
get<
bool>(
"aggregation: brick x Dirichlet");
167 dirichletY_ = pL.
get<
bool>(
"aggregation: brick y Dirichlet");
168 dirichletZ_ = pL.
get<
bool>(
"aggregation: brick z Dirichlet");
169 if (dirichletX_) GetOStream(
Runtime0) <<
"Dirichlet boundaries in the x direction" << std::endl;
170 if (dirichletY_) GetOStream(
Runtime0) <<
"Dirichlet boundaries in the y direction" << std::endl;
171 if (dirichletZ_) GetOStream(
Runtime0) <<
"Dirichlet boundaries in the z direction" << std::endl;
187 Setup(comm, overlappedCoords, colMap);
189 GetOStream(
Runtime0) <<
"Using brick size: " << bx_
190 << (nDim_ > 1 ?
"x " +
toString(by_) :
"")
191 << (nDim_ > 2 ?
"x " +
toString(bz_) :
"") << std::endl;
194 BuildGraph(currentLevel, A);
210 std::set<GO> myAggGIDs, remoteAggGIDs;
211 for (LO LID = 0; LID < numPoints; LID++) {
212 GO aggGID = getAggGID(LID);
214 if (aggGID == GO_INVALID)
continue;
217 if ((revMap_.find(getRoot(LID)) != revMap_.end()) && rowMap->isNodeGlobalElement(colMap->getGlobalElement(revMap_[getRoot(LID)]))) {
219 vertex2AggId[LID] = aggGID;
220 myAggGIDs.insert(aggGID);
226 remoteAggGIDs.insert(aggGID);
229 size_t numAggregates = myAggGIDs.size();
230 size_t numRemote = remoteAggGIDs.size();
233 std::map<GO, LO> AggG2L;
234 std::map<GO, int> AggG2R;
236 Array<GO> myAggGIDsArray(numAggregates), remoteAggGIDsArray(numRemote);
240 for (
typename std::set<GO>::const_iterator it = myAggGIDs.begin(); it != myAggGIDs.end(); it++) {
242 AggG2R[*it] = myRank;
244 myAggGIDsArray[ind++] = *it;
249 myAggGIDsArray, 0, comm);
252 for (
typename std::set<GO>::const_iterator it = remoteAggGIDs.begin(); it != remoteAggGIDs.end(); it++)
253 remoteAggGIDsArray[ind++] = *it;
258 aggMap->getRemoteIndexList(remoteAggGIDsArray, remoteProcIDs, remoteLIDs);
261 for (
size_t i = 0; i < numRemote; i++) {
262 AggG2L[remoteAggGIDsArray[i]] = remoteLIDs[i];
263 AggG2R[remoteAggGIDsArray[i]] = remoteProcIDs[i];
267 for (LO LID = 0; LID < numPoints; LID++) {
268 if (revMap_.find(getRoot(LID)) != revMap_.end() && rowMap->isNodeGlobalElement(colMap->getGlobalElement(revMap_[getRoot(LID)]))) {
269 GO aggGID = vertex2AggId[LID];
271 vertex2AggId[LID] = AggG2L[aggGID];
272 procWinner[LID] = AggG2R[aggGID];
281 Set(currentLevel,
"Aggregates", aggregates);
286 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
289 nDim_ = coords->getNumVectors();
291 x_ = coords->getData(0);
292 xMap_ = Construct1DMap(comm, x_);
297 y_ = coords->getData(1);
298 yMap_ = Construct1DMap(comm, y_);
304 z_ = coords->getData(2);
305 zMap_ = Construct1DMap(comm, z_);
309 for (
size_t ind = 0; ind < coords->getLocalLength(); ind++) {
310 GO i = (*xMap_)[(coords->getData(0))[ind]], j = 0, k = 0;
312 j = (*yMap_)[(coords->getData(1))[ind]];
314 k = (*zMap_)[(coords->getData(2))[ind]];
316 revMap_[k * ny_ * nx_ + j * nx_ + i] = ind;
320 int xboost = dirichletX_ ? 1 : 0;
321 int yboost = dirichletY_ ? 1 : 0;
322 int zboost = dirichletZ_ ? 1 : 0;
323 naggx_ = (nx_ - 2 * xboost) / bx_ + ((nx_ - 2 * xboost) % bx_ ? 1 : 0);
326 naggy_ = (ny_ - 2 * yboost) / by_ + ((ny_ - 2 * yboost) % by_ ? 1 : 0);
331 naggz_ = (nz_ - 2 * zboost) / bz_ + ((nz_ - 2 * zboost) % bz_ ? 1 : 0);
336 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
345 for (
int i = 0; i < n; i++)
352 int numProcs = comm->getSize();
356 MPI_Comm rawComm = (*dupMpiComm->getRawMpiComm())();
358 int sendCnt = gMap->size(), cnt = 0, recvSize;
359 Array<int> recvCnt(numProcs), Displs(numProcs);
363 for (
typename container::const_iterator cit = gMap->begin(); cit != gMap->end(); cit++)
364 sendBuf[cnt++] = Teuchos::as<double>(STS::real(cit->first));
366 MPI_Allgather(&sendCnt, 1, MPI_INT, recvCnt.getRawPtr(), 1, MPI_INT, rawComm);
368 for (
int i = 0; i < numProcs - 1; i++)
369 Displs[i + 1] = Displs[i] + recvCnt[i];
370 recvSize = Displs[numProcs - 1] + recvCnt[numProcs - 1];
372 MPI_Allgatherv(sendBuf.
getRawPtr(), sendCnt, MPI_DOUBLE, recvBuf.
getRawPtr(), recvCnt.getRawPtr(), Displs.getRawPtr(), MPI_DOUBLE, rawComm);
374 for (
int i = 0; i < recvSize; i++)
375 (*gMap)[as<SC>(recvBuf[i])] = 0;
380 for (
typename container::iterator it = gMap->begin(); it != gMap->end(); it++)
386 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
389 getIJK(LID, i, j, k);
391 return (k * ny_ * nx_ + j * nx_ + i) == getRoot(LID);
394 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
396 bool boundary =
false;
398 getIJK(LID, i, j, k);
399 if (dirichletX_ && (i == 0 || i == nx_ - 1))
401 if (nDim_ > 1 && dirichletY_ && (j == 0 || j == ny_ - 1))
403 if (nDim_ > 2 && dirichletZ_ && (k == 0 || k == nz_ - 1))
409 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
411 if (isDirichlet(LID))
414 int aggI, aggJ, aggK;
415 getAggIJK(LID, aggI, aggJ, aggK);
416 int xboost = dirichletX_ ? 1 : 0;
417 int yboost = dirichletY_ ? 1 : 0;
418 int zboost = dirichletZ_ ? 1 : 0;
420 int i = xboost + aggI * bx_ + (bx_ - 1) / 2;
421 int j = (nDim_ > 1) ? yboost + aggJ * by_ + (by_ - 1) / 2 : 0;
422 int k = (nDim_ > 2) ? zboost + aggK * bz_ + (bz_ - 1) / 2 : 0;
424 return k * ny_ * nx_ + j * nx_ + i;
427 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
429 i = (*xMap_)[x_[LID]];
430 j = (nDim_ > 1) ? (*yMap_)[y_[LID]] : 0;
431 k = (nDim_ > 2) ? (*zMap_)[z_[LID]] : 0;
434 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
436 int xboost = dirichletX_ ? 1 : 0;
437 int yboost = dirichletY_ ? 1 : 0;
438 int zboost = dirichletZ_ ? 1 : 0;
439 int pointI, pointJ, pointK;
440 getIJK(LID, pointI, pointJ, pointK);
441 i = (pointI - xboost) / bx_;
444 j = (pointJ - yboost) / by_;
449 k = (pointK - zboost) / bz_;
454 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
456 bool boundary =
false;
459 getIJK(LID, i, j, k);
461 getAggIJK(LID, ii, jj, kk);
463 if (dirichletX_ && (i == 0 || i == nx_ - 1)) boundary =
true;
464 if (nDim_ > 1 && dirichletY_ && (j == 0 || j == ny_ - 1)) boundary =
true;
465 if (nDim_ > 2 && dirichletZ_ && (k == 0 || k == nz_ - 1)) boundary =
true;
477 return Teuchos::as<GlobalOrdinal>(kk * naggy_ * naggx_) + Teuchos::as<GlobalOrdinal>(jj * naggx_) + ii;
480 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
483 double dirichletThreshold = 0.0;
485 if (bx_ > 1 && (nDim_ <= 1 || by_ > 1) && (nDim_ <= 2 || bz_ > 1)) {
486 FactoryMonitor m(*
this,
"Generating Graph (trivial)", currentLevel);
494 GO numLocalBoundaryNodes = 0;
495 GO numGlobalBoundaryNodes = 0;
496 for (
size_t i = 0; i < boundaryNodes.size(); ++i)
497 if (boundaryNodes(i))
498 numLocalBoundaryNodes++;
500 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
501 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes <<
" Dirichlet nodes" << std::endl;
503 Set(currentLevel,
"DofsPerNode", 1);
504 Set(currentLevel,
"Graph", graph);
505 Set(currentLevel,
"Filtering",
false);
511 bool drop_x = (bx_ == 1);
512 bool drop_y = (nDim_ > 1 && by_ == 1);
513 bool drop_z = (nDim_ > 2 && bz_ == 1);
515 typename LWGraph::row_type::non_const_type
rows(
"rows", A->getLocalNumRows() + 1);
516 typename LWGraph::entries_type::non_const_type columns(
"columns", A->getLocalNumEntries());
518 size_t N = A->getRowMap()->getLocalNumElements();
521 auto G = A->getLocalMatrixHost().graph;
522 auto rowptr = G.row_map;
523 auto colind = G.entries;
527 for (
size_t row = 0; row < N; row++) {
530 LO row2 = A->getColMap()->getLocalElement(A->getRowMap()->getGlobalElement(row));
531 getIJK(row2, ir, jr, kr);
533 for (
size_t cidx = rowptr[row]; cidx < rowptr[row + 1]; cidx++) {
535 LO col = colind[cidx];
536 getIJK(col, ic, jc, kc);
538 if ((row2 != col) && ((drop_x && ir != ic) || (drop_y && jr != jc) || (drop_z && kr != kc))) {
557 GO numLocalBoundaryNodes = 0;
558 GO numGlobalBoundaryNodes = 0;
559 for (
size_t i = 0; i < boundaryNodes.size(); ++i)
560 if (boundaryNodes(i))
561 numLocalBoundaryNodes++;
563 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
564 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes <<
" Dirichlet nodes" << std::endl;
566 Set(currentLevel,
"DofsPerNode", 1);
567 Set(currentLevel,
"Graph", graph);
568 Set(currentLevel,
"Filtering",
true);
#define MueLu_sumAll(rcpComm, in, out)
void getAggIJK(LocalOrdinal LID, int &i, int &j, int &k) const
MueLu::DefaultLocalOrdinal LocalOrdinal
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
const RCP< LOVector > & GetProcWinner() const
Returns constant vector that maps local node IDs to owning processor IDs.
KOKKOS_INLINE_FUNCTION void SetBoundaryNodeMap(const boundary_nodes_type bndry)
Set boolean array indicating which rows correspond to Dirichlet boundaries.
Container class for aggregation information.
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.
One-liner description of what is happening.
void Build(Level ¤tLevel) const
Build aggregates.
void SetIsRoot(LO i, bool value=true)
Set root node information.
GlobalOrdinal getAggGID(LocalOrdinal LID) const
std::map< Scalar, GlobalOrdinal, compare > container
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Class that holds all level-specific information.
void Setup(const RCP< const Teuchos::Comm< int > > &comm, const RCP< Xpetra::MultiVector< typename Teuchos::ScalarTraits< Scalar >::magnitudeType, LO, GO, NO > > &coords, const RCP< const Map > &map) const
bool isDirichlet(LocalOrdinal LID) const
#define MUELU_UNAGGREGATED
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 resize(size_type new_size, const value_type &x=value_type())
bool isRoot(LocalOrdinal LID) const
void getIJK(LocalOrdinal LID, int &i, int &j, int &k) const
Lightweight MueLu representation of a compressed row storage graph.
RCP< container > Construct1DMap(const RCP< const Teuchos::Comm< int > > &comm, const ArrayRCP< const typename Teuchos::ScalarTraits< Scalar >::magnitudeType > &x) const
KOKKOS_INLINE_FUNCTION void AggregatesCrossProcessors(const bool &flag)
Record whether aggregates include DOFs from other processes.
static Kokkos::View< bool *, typename Kokkos::HostSpace > DetectDirichletRows_kokkos_host(const Matrix &A, const Magnitude &tol=Teuchos::ScalarTraits< typename Teuchos::ScalarTraits< SC >::magnitudeType >::zero(), const bool count_twos_as_dirichlet=false)
std::string description() const
Return a simple one-line description of this object.
void BuildGraph(Level ¤tLevel, const RCP< Matrix > &A) const
GlobalOrdinal getRoot(LocalOrdinal LID) const
void DeclareInput(Level ¤tLevel) const
Input.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
void SetNumAggregates(LO nAggregates)
Set number of local aggregates on current processor.
#define SET_VALID_ENTRY(name)