46 #ifndef MUELU_COORDINATESTRANSFER_FACTORY_DEF_HPP 47 #define MUELU_COORDINATESTRANSFER_FACTORY_DEF_HPP 49 #include "Xpetra_ImportFactory.hpp" 50 #include "Xpetra_MultiVectorFactory.hpp" 51 #include "Xpetra_MapFactory.hpp" 52 #include "Xpetra_IO.hpp" 54 #include "MueLu_CoarseMapFactory.hpp" 55 #include "MueLu_Aggregates.hpp" 64 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
66 RCP<ParameterList> validParamList = rcp(
new ParameterList());
68 validParamList->set< RCP<const FactoryBase> >(
"Coordinates", Teuchos::null,
"Factory for coordinates generation");
69 validParamList->set< RCP<const FactoryBase> >(
"Aggregates", Teuchos::null,
"Factory for coordinates generation");
70 validParamList->set< RCP<const FactoryBase> >(
"CoarseMap", Teuchos::null,
"Generating factory of the coarse map");
71 validParamList->set<
int > (
"write start", -1,
"first level at which coordinates should be written to file");
72 validParamList->set<
int > (
"write end", -1,
"last level at which coordinates should be written to file");
74 return validParamList;
77 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
79 static bool isAvailableCoords =
false;
82 isAvailableCoords = coarseLevel.
IsAvailable(
"Coordinates",
this);
84 if (isAvailableCoords ==
false) {
85 Input(fineLevel,
"Coordinates");
86 Input(fineLevel,
"Aggregates");
87 Input(fineLevel,
"CoarseMap");
91 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
95 typedef Xpetra::MultiVector<double,LO,GO,NO> xdMV;
97 GetOStream(
Runtime0) <<
"Transferring coordinates" << std::endl;
100 GetOStream(
Runtime0) <<
"Reusing coordinates" << std::endl;
104 RCP<Aggregates> aggregates = Get< RCP<Aggregates> > (fineLevel,
"Aggregates");
105 RCP<xdMV> fineCoords = Get< RCP<xdMV> >(fineLevel,
"Coordinates");
106 RCP<const Map> coarseMap = Get< RCP<const Map> > (fineLevel,
"CoarseMap");
112 ArrayView<const GO> elementAList = coarseMap->getNodeElementList();
114 if (rcp_dynamic_cast<const StridedMap>(coarseMap) != Teuchos::null)
115 blkSize = rcp_dynamic_cast<
const StridedMap>(coarseMap)->getFixedBlockSize();
117 GO indexBase = coarseMap->getIndexBase();
118 size_t numElements = elementAList.size() / blkSize;
119 Array<GO> elementList(numElements);
122 for (LO i = 0; i < Teuchos::as<LO>(numElements); i++)
123 elementList[i] = (elementAList[i*blkSize]-indexBase)/blkSize + indexBase;
125 RCP<const Map> uniqueMap = fineCoords->getMap();
126 RCP<const Map> coarseCoordMap = MapFactory ::Build(coarseMap->lib(), Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(), elementList, indexBase, coarseMap->getComm());
127 RCP<xdMV> coarseCoords = Xpetra::MultiVectorFactory<double,LO,GO,NO>::Build(coarseCoordMap, fineCoords->getNumVectors());
130 RCP<xdMV> ghostedCoords = fineCoords;
131 if (aggregates->AggregatesCrossProcessors()) {
132 RCP<const Map> nonUniqueMap = aggregates->GetMap();
133 RCP<const Import> importer = ImportFactory::Build(uniqueMap, nonUniqueMap);
135 ghostedCoords = Xpetra::MultiVectorFactory<double,LO,GO,NO>::Build(nonUniqueMap, fineCoords->getNumVectors());
136 ghostedCoords->doImport(*fineCoords, *importer, Xpetra::INSERT);
140 int myPID = uniqueMap->getComm()->getRank();
141 LO numAggs = aggregates->GetNumAggregates();
142 ArrayRCP<LO> aggSizes = aggregates->ComputeAggregateSizes(
true,
true);
143 const ArrayRCP<const LO> vertex2AggID = aggregates->GetVertex2AggId()->getData(0);
144 const ArrayRCP<const LO> procWinner = aggregates->GetProcWinner()->getData(0);
147 for (
size_t j = 0; j < fineCoords->getNumVectors(); j++) {
148 ArrayRCP<const double> fineCoordsData = ghostedCoords->getData(j);
149 ArrayRCP<double> coarseCoordsData = coarseCoords->getDataNonConst(j);
151 for (LO lnode = 0; lnode < vertex2AggID.size(); lnode++)
152 if (procWinner[lnode] == myPID)
153 coarseCoordsData[vertex2AggID[lnode]] += fineCoordsData[lnode];
155 for (LO agg = 0; agg < numAggs; agg++)
156 coarseCoordsData[agg] /= aggSizes[agg];
159 Set<RCP<xdMV> >(coarseLevel,
"Coordinates", coarseCoords);
161 const ParameterList& pL = GetParameterList();
162 int writeStart = pL.get<
int>(
"write start"), writeEnd = pL.get<
int>(
"write end");
163 if (writeStart == 0 && fineLevel.GetLevelID() == 0 && writeStart <= writeEnd) {
164 std::ostringstream buf;
165 buf << fineLevel.GetLevelID();
166 std::string fileName =
"coordinates_before_rebalance_level_" + buf.str() +
".m";
167 Xpetra::IO<double,LO,GO,NO>::Write(fileName,*fineCoords);
169 if (writeStart <= coarseLevel.GetLevelID() && coarseLevel.GetLevelID() <= writeEnd) {
170 std::ostringstream buf;
171 buf << coarseLevel.GetLevelID();
172 std::string fileName =
"coordinates_before_rebalance_level_" + buf.str() +
".m";
173 Xpetra::IO<double,LO,GO,NO>::Write(fileName,*coarseCoords);
179 #endif // MUELU_COORDINATESTRANSFER_FACTORY_DEF_HPP void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
void DeclareInput(Level &finelevel, Level &coarseLevel) const
Specifies the data that this class needs, and the factories that generate that data.
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.
RequestMode GetRequestMode() const
Class that holds all level-specific information.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.