46 #ifndef MUELU_BRICKAGGREGATIONFACTORY_DEF_HPP_ 47 #define MUELU_BRICKAGGREGATIONFACTORY_DEF_HPP_ 51 #include <Teuchos_DefaultMpiComm.hpp> 52 #include <Teuchos_CommHelpers.hpp> 54 #include <Teuchos_OrdinalTraits.hpp> 56 #include <Xpetra_Import.hpp> 57 #include <Xpetra_ImportFactory.hpp> 58 #include <Xpetra_Map.hpp> 59 #include <Xpetra_MapFactory.hpp> 60 #include <Xpetra_Matrix.hpp> 61 #include <Xpetra_MultiVector.hpp> 62 #include <Xpetra_MultiVectorFactory.hpp> 65 #include "MueLu_Aggregates.hpp" 69 #include "MueLu_Utilities.hpp" 74 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
76 RCP<ParameterList> validParamList = rcp(
new ParameterList());
78 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name)) 85 #undef SET_VALID_ENTRY 87 validParamList->set< RCP<const FactoryBase> >(
"A", Teuchos::null,
"Generating factory for matrix");
88 validParamList->set< RCP<const FactoryBase> >(
"Coordinates", Teuchos::null,
"Generating factory for coordinates");
89 validParamList->set< RCP<const FactoryBase> >(
"Graph" , Teuchos::null,
"Generating factory of the graph");
90 validParamList->set< RCP<const FactoryBase> >(
"DofsPerNode", null,
"Generating factory for variable \'DofsPerNode\', usually the same as for \'Graph\'");
92 return validParamList;
95 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
97 Input(currentLevel,
"A");
98 Input(currentLevel,
"Coordinates");
145 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
149 typedef Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO> MultiVector_d;
151 const ParameterList& pL = GetParameterList();
152 RCP<MultiVector_d> coords = Get<RCP<MultiVector_d> >(currentLevel,
"Coordinates");
153 RCP<Matrix> A = Get< RCP<Matrix> > (currentLevel,
"A");
154 RCP<const Map> rowMap = A->getRowMap();
155 RCP<const Map> colMap = A->getColMap();
156 GO GO_INVALID = Teuchos::OrdinalTraits<GO>::invalid();
158 RCP<const Teuchos::Comm<int> > comm = rowMap->getComm();
159 int numProcs = comm->getSize();
160 int myRank = comm->getRank();
162 int numPoints = colMap->getNodeNumElements();
164 bx_ = pL.get<
int>(
"aggregation: brick x size");
165 by_ = pL.get<
int>(
"aggregation: brick y size");
166 bz_ = pL.get<
int>(
"aggregation: brick z size");
168 dirichletX_ = pL.get<
bool>(
"aggregation: brick x Dirichlet");
169 dirichletY_ = pL.get<
bool>(
"aggregation: brick y Dirichlet");
170 dirichletZ_ = pL.get<
bool>(
"aggregation: brick z Dirichlet");
171 if(dirichletX_) GetOStream(
Runtime0) <<
"Dirichlet boundaries in the x direction"<<std::endl;
172 if(dirichletY_) GetOStream(
Runtime0) <<
"Dirichlet boundaries in the y direction"<<std::endl;
173 if(dirichletZ_) GetOStream(
Runtime0) <<
"Dirichlet boundaries in the z direction"<<std::endl;
180 RCP<MultiVector_d> overlappedCoords = coords;
181 RCP<const Import> importer = ImportFactory::Build(coords->getMap(), colMap);
182 if (!importer.is_null()) {
183 overlappedCoords = Xpetra::MultiVectorFactory<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO>::Build(colMap, coords->getNumVectors());
184 overlappedCoords->doImport(*coords, *importer, Xpetra::INSERT);
189 Setup(comm, overlappedCoords, colMap);
191 GetOStream(
Runtime0) <<
"Using brick size: " << bx_
192 << (nDim_ > 1 ?
"x " +
toString(by_) :
"")
193 << (nDim_ > 2 ?
"x " +
toString(bz_) :
"") << std::endl;
196 RCP<Aggregates> aggregates = rcp(
new Aggregates(colMap));
197 aggregates->setObjectLabel(
"Brick");
199 ArrayRCP<LO> vertex2AggId = aggregates->GetVertex2AggId()->getDataNonConst(0);
200 ArrayRCP<LO> procWinner = aggregates->GetProcWinner() ->getDataNonConst(0);
209 std::set<GO> myAggGIDs, remoteAggGIDs;
210 for (LO LID = 0; LID < numPoints; LID++) {
211 GO aggGID = getAggGID(LID);
213 if(aggGID == GO_INVALID)
continue;
216 if ((revMap_.find(getRoot(LID)) != revMap_.end()) && rowMap->isNodeGlobalElement(colMap->getGlobalElement(revMap_[getRoot(LID)]))) {
218 vertex2AggId[LID] = aggGID;
219 myAggGIDs.insert(aggGID);
222 aggregates->SetIsRoot(LID);
225 remoteAggGIDs.insert(aggGID);
228 size_t numAggregates = myAggGIDs .size();
229 size_t numRemote = remoteAggGIDs.size();
230 aggregates->SetNumAggregates(numAggregates);
232 std::map<GO,LO> AggG2L;
233 std::map<GO,int> AggG2R;
235 Array<GO> myAggGIDsArray(numAggregates), remoteAggGIDsArray(numRemote);
239 for (
typename std::set<GO>::const_iterator it = myAggGIDs.begin(); it != myAggGIDs.end(); it++) {
241 AggG2R[*it] = myRank;
243 myAggGIDsArray[ind++] = *it;
247 RCP<Map> aggMap = MapFactory::Build(rowMap->lib(), Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
248 myAggGIDsArray, 0, comm);
251 for (
typename std::set<GO>::const_iterator it = remoteAggGIDs.begin(); it != remoteAggGIDs.end(); it++)
252 remoteAggGIDsArray[ind++] = *it;
255 Array<int> remoteProcIDs(numRemote);
256 Array<LO> remoteLIDs (numRemote);
257 aggMap->getRemoteIndexList(remoteAggGIDsArray, remoteProcIDs, remoteLIDs);
260 for (
size_t i = 0; i < numRemote; i++) {
261 AggG2L[remoteAggGIDsArray[i]] = remoteLIDs [i];
262 AggG2R[remoteAggGIDsArray[i]] = remoteProcIDs[i];
266 for (LO LID = 0; LID < numPoints; LID++) {
267 if (revMap_.find(getRoot(LID)) != revMap_.end() && rowMap->isNodeGlobalElement(colMap->getGlobalElement(revMap_[getRoot(LID)]))) {
268 GO aggGID = vertex2AggId[LID];
270 vertex2AggId[LID] = AggG2L[aggGID];
271 procWinner [LID] = AggG2R[aggGID];
279 aggregates->AggregatesCrossProcessors(numGlobalRemote);
281 Set(currentLevel,
"Aggregates", aggregates);
283 GetOStream(
Statistics1) << aggregates->description() << std::endl;
286 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
288 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>& )
const {
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;
321 int xboost = dirichletX_ ? 1 : 0;
322 int yboost = dirichletY_ ? 1 : 0;
323 int zboost = dirichletZ_ ? 1 : 0;
324 naggx_ = (nx_-2*xboost)/bx_ + ((nx_-2*xboost) % bx_ ? 1 : 0);
327 naggy_ = (ny_-2*yboost)/by_ + ( (ny_-2*yboost) % by_ ? 1 : 0);
332 naggz_ = (nz_-2*zboost)/bz_ + ( (nz_-2*zboost) % bz_ ? 1 : 0);
338 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
339 RCP<typename BrickAggregationFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::container>
342 const ArrayRCP<
const typename Teuchos::ScalarTraits<Scalar>::magnitudeType>& x)
const 347 RCP<container> gMap = rcp(
new container);
348 for (
int i = 0; i < n; i++)
355 int numProcs = comm->getSize();
357 RCP<const Teuchos::MpiComm<int> > dupMpiComm = rcp_dynamic_cast<
const Teuchos::MpiComm<int> >(comm->duplicate());
359 MPI_Comm rawComm = (*dupMpiComm->getRawMpiComm())();
361 int sendCnt = gMap->size(), cnt = 0, recvSize;
362 Array<int> recvCnt(numProcs), Displs(numProcs);
363 Array<double> sendBuf, recvBuf;
365 sendBuf.resize(sendCnt);
366 for (
typename container::const_iterator cit = gMap->begin(); cit != gMap->end(); cit++)
367 sendBuf[cnt++] = Teuchos::as<double>(STS::real(cit->first));
369 MPI_Allgather(&sendCnt, 1, MPI_INT, recvCnt.getRawPtr(), 1, MPI_INT, rawComm);
371 for (
int i = 0; i < numProcs-1; i++)
372 Displs[i+1] = Displs[i] + recvCnt[i];
373 recvSize = Displs[numProcs-1] + recvCnt[numProcs-1];
374 recvBuf.resize(recvSize);
375 MPI_Allgatherv(sendBuf.getRawPtr(), sendCnt, MPI_DOUBLE, recvBuf.getRawPtr(), recvCnt.getRawPtr(), Displs.getRawPtr(), MPI_DOUBLE, rawComm);
377 for (
int i = 0; i < recvSize; i++)
378 (*gMap)[as<SC>(recvBuf[i])] = 0;
383 for (
typename container::iterator it = gMap->begin(); it != gMap->end(); it++)
389 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
394 return (k*ny_*nx_ + j*nx_ + i) == getRoot(LID);
397 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
399 bool boundary =
false;
402 if( dirichletX_ && (i == 0 || i == nx_-1) )
404 if(nDim_ > 1 && dirichletY_ && (j == 0 || j == ny_-1) )
406 if(nDim_ > 2 && dirichletZ_ && (k == 0 || k == nz_-1) )
413 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
416 return Teuchos::OrdinalTraits<GlobalOrdinal>::invalid();
419 getAggIJK(LID,aggI,aggJ,aggK);
420 int xboost = dirichletX_ ? 1 : 0;
421 int yboost = dirichletY_ ? 1 : 0;
422 int zboost = dirichletZ_ ? 1 : 0;
424 int i = xboost + aggI*bx_ + (bx_-1)/2;
425 int j = (nDim_>1) ? yboost + aggJ*by_ + (by_-1)/2 : 0;
426 int k = (nDim_>2) ? zboost + aggK*bz_ + (bz_-1)/2 : 0;
428 return k*ny_*nx_ + j*nx_ + i;
431 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
433 i = (*xMap_)[x_[LID]];
434 j = (nDim_>1) ? (*yMap_)[y_[LID]] : 0;
435 k = (nDim_>2) ? (*zMap_)[z_[LID]] : 0;
439 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
441 int xboost = dirichletX_ ? 1 : 0;
442 int yboost = dirichletY_ ? 1 : 0;
443 int zboost = dirichletZ_ ? 1 : 0;
444 int pointI, pointJ, pointK;
445 getIJK(LID,pointI,pointJ,pointK);
446 i = (pointI-xboost)/bx_;
448 if (nDim_ > 1) j = (pointJ-yboost)/by_;
451 if (nDim_ > 2) k = (pointK-zboost)/bz_;
455 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
457 bool boundary =
false;
462 getAggIJK(LID,ii,jj,kk);
464 if( dirichletX_ && (i == 0 || i == nx_ - 1)) boundary =
true;
465 if (nDim_ > 1 && dirichletY_ && (j == 0 || j == ny_ - 1)) boundary =
true;
466 if (nDim_ > 2 && dirichletZ_ && (k == 0 || k == nz_ - 1)) boundary =
true;
476 return Teuchos::OrdinalTraits<GlobalOrdinal>::invalid();
478 return Teuchos::as<GlobalOrdinal>(kk*naggy_*naggx_) + Teuchos::as<GlobalOrdinal>(jj*naggx_) + ii;
#define MueLu_sumAll(rcpComm, in, out)
MueLu::DefaultLocalOrdinal LocalOrdinal
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
Container class for aggregation information.
void DeclareInput(Level ¤tLevel) const
Input.
void getAggIJK(LocalOrdinal LID, int &i, int &j, int &k) const
Timer to be used in factories. Similar to Monitor but with additional timers.
One-liner description of what is happening.
Namespace for MueLu classes and methods.
bool isDirichlet(LocalOrdinal LID) const
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Class that holds all level-specific information.
GlobalOrdinal getRoot(LocalOrdinal LID) const
#define MUELU_UNAGGREGATED
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
RCP< container > Construct1DMap(const RCP< const Teuchos::Comm< int > > &comm, const ArrayRCP< const typename Teuchos::ScalarTraits< Scalar >::magnitudeType > &x) const
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
void Build(Level ¤tLevel) const
Build aggregates.
bool isRoot(LocalOrdinal LID) const
std::map< Scalar, GlobalOrdinal, compare > container
void getIJK(LocalOrdinal LID, int &i, int &j, int &k) const
#define SET_VALID_ENTRY(name)
GlobalOrdinal getAggGID(LocalOrdinal LID) const