47 #ifndef MUELU_COARSENINGVISUALIZATIONFACTORY_DEF_HPP_
48 #define MUELU_COARSENINGVISUALIZATIONFACTORY_DEF_HPP_
52 #include <Xpetra_MultiVectorFactory.hpp>
58 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
62 validParamList->
set<
int>(
"visualization: start level", 0,
"visualize only levels with level ids greater or equal than start level");
64 validParamList->
set<
RCP<const FactoryBase> >(
"P", Teuchos::null,
"Prolongator factory. The user has to declare either P or Ptent but not both at the same time.");
65 validParamList->
set<
RCP<const FactoryBase> >(
"Ptent", Teuchos::null,
"Tentative prolongator factory. The user has to declare either P or Ptent as input but not both at the same time");
69 return validParamList;
72 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
74 this->Input(fineLevel,
"Coordinates");
78 pL.
isParameter(
"Ptent") && GetFactory(
"Ptent") != Teuchos::null,
80 "You must not declare both P and Ptent. Use only once for visualization.");
82 "You have to either declare P or Ptent for visualization, but not both.");
84 if (GetFactory(
"P") != Teuchos::null && GetFactory(
"Ptent") == Teuchos::null)
85 this->Input(coarseLevel,
"P");
86 else if (GetFactory(
"Ptent") != Teuchos::null && GetFactory(
"P") == Teuchos::null)
87 this->Input(coarseLevel,
"Ptent");
89 if (pL.
get<
bool>(
"visualization: fine graph edges"))
90 Input(fineLevel,
"Graph");
92 if(pL.
get<
bool>(
"visualization: coarse graph edges")) {
93 Input(coarseLevel,
"Coordinates");
94 Input(coarseLevel,
"Graph");
99 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
104 if (this->GetFactory(
"P") != Teuchos::null && this->GetFactory(
"Ptent") == Teuchos::null)
105 P = Get<RCP<Matrix> >(coarseLevel,
"P");
106 if (GetFactory(
"Ptent") != Teuchos::null && GetFactory(
"P") == Teuchos::null)
114 if (P->IsView(
"stridedMaps") && Teuchos::rcp_dynamic_cast<
const StridedMap>(P->getRowMap(
"stridedMaps")) != Teuchos::null) {
115 strRowMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(P->getRowMap(
"stridedMaps"));
118 std::vector<size_t> stridingInfo = strRowMap->getStridingData();
119 for (
size_t j = 0; j < Teuchos::as<size_t>(blockid); j++)
120 stridedRowOffset += stridingInfo[j];
121 dofsPerNode = Teuchos::as<LocalOrdinal>(stridingInfo[blockid]);
123 dofsPerNode = strRowMap->getFixedBlockSize();
125 GetOStream(
Runtime1) <<
"CoarseningVisualizationFactory::Build():"
126 <<
" #dofs per node = " << dofsPerNode << std::endl;
132 if (P->IsView(
"stridedMaps") && Teuchos::rcp_dynamic_cast<
const StridedMap>(P->getRowMap(
"stridedMaps")) != Teuchos::null) {
133 strDomainMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(P->getColMap(
"stridedMaps"));
134 LocalOrdinal blockid = strDomainMap->getStridedBlockId();
137 std::vector<size_t> stridingInfo = strDomainMap->getStridingData();
138 for (
size_t j = 0; j < Teuchos::as<size_t>(blockid); j++)
139 stridedColumnOffset += stridingInfo[j];
140 columnsPerNode = Teuchos::as<LocalOrdinal>(stridingInfo[blockid]);
142 columnsPerNode = strDomainMap->getFixedBlockSize();
144 GetOStream(
Runtime1) <<
"CoarseningVisualizationFactory::Build():"
145 <<
" #columns per node = " << columnsPerNode << std::endl;
150 "CoarseningVisualization only supports non-overlapping transfers");
153 LocalOrdinal numLocalAggs = strDomainMap->getLocalNumElements() / columnsPerNode;
154 std::vector<std::set<LocalOrdinal> > localAggs(numLocalAggs);
157 for (
LO row = 0; row < Teuchos::as<LO>(P->getRowMap()->getLocalNumElements()); ++row) {
160 P->getLocalRowView(row, indices, vals);
163 localAggs[(*c) / columnsPerNode].insert(row / dofsPerNode);
168 std::vector<int> myLocalAggsPerProc(comm->getSize(), 0);
169 std::vector<int> numLocalAggsPerProc(comm->getSize(), 0);
170 myLocalAggsPerProc[comm->getRank()] = numLocalAggs;
171 Teuchos::reduceAll(*comm,
Teuchos::REDUCE_MAX, comm->getSize(), &myLocalAggsPerProc[0], &numLocalAggsPerProc[0]);
174 for (
int i = 0; i < comm->getRank(); ++i) {
175 myAggOffset += numLocalAggsPerProc[i];
191 "Number of fine level nodes in coordinates is inconsistent with dof based information");
194 if (pL.
get<
bool>(
"visualization: fine graph edges")) {
195 fineGraph = Get<RCP<LWGraph> >(fineLevel,
"Graph");
200 coords = ghostedCoords;
208 if (coords->getNumVectors() == 3) {
213 LocalOrdinal numFineNodes = Teuchos::as<LocalOrdinal>(coords->getLocalLength());
217 for (LocalOrdinal i = 0; i < numLocalAggs; ++i) {
219 for (
typename std::set<LocalOrdinal>::iterator it = localAggs[i].begin(); it != localAggs[i].end(); ++it) {
220 vertex2AggId[*it] = i;
227 std::vector<bool> isRoot(numFineNodes,
false);
228 for (LocalOrdinal i = 0; i < numLocalAggs; ++i) {
234 for (
typename std::set<LocalOrdinal>::iterator it = localAggs[i].begin(); it != localAggs[i].end(); ++it) {
235 xCenter += xCoords[*it];
236 yCenter += yCoords[*it];
237 if (coords->getNumVectors() == 3) zCenter += zCoords[*it];
239 xCenter /= Teuchos::as<LocalOrdinal>(localAggs[i].size());
240 yCenter /= Teuchos::as<LocalOrdinal>(localAggs[i].size());
241 zCenter /= Teuchos::as<LocalOrdinal>(localAggs[i].size());
244 LocalOrdinal rootCandidate = -1;
246 for (
typename std::set<LocalOrdinal>::iterator it = localAggs[i].begin(); it != localAggs[i].end(); ++it) {
250 if (coords->getNumVectors() == 3) tempz = zCenter - zCoords[*it];
252 mydistance += tempx * tempx;
253 mydistance += tempy * tempy;
254 mydistance += tempz * tempz;
255 mydistance = sqrt(mydistance);
256 if (mydistance <= minDistance) {
257 minDistance = mydistance;
262 isRoot[rootCandidate] =
true;
265 std::vector<LocalOrdinal> vertices;
266 std::vector<LocalOrdinal> geomSize;
267 int vizLevel = pL.
get<
int>(
"visualization: start level");
268 if (vizLevel <= fineLevel.GetLevelID()) {
269 std::string aggStyle = pL.
get<std::string>(
"visualization: style");
270 if (aggStyle ==
"Point Cloud")
271 this->doPointCloud(vertices, geomSize, numLocalAggs, numFineNodes);
272 else if (aggStyle ==
"Jacks")
273 this->doJacks(vertices, geomSize, numLocalAggs, numFineNodes, isRoot, vertex2AggId);
274 else if (aggStyle ==
"Convex Hulls") {
277 if (coords->getNumVectors() == 3)
278 this->doConvexHulls3D(vertices, geomSize, numLocalAggs, numFineNodes, isRoot, vertex2AggId, xCoords, yCoords, zCoords);
279 else if (coords->getNumVectors() == 2)
280 this->doConvexHulls2D(vertices, geomSize, numLocalAggs, numFineNodes, isRoot, vertex2AggId, xCoords, yCoords);
282 GetOStream(
Warnings0) <<
" Warning: Unrecognized agg style.\nPossible values are Point Cloud, Jacks, Convex Hulls.\nDefaulting to Point Cloud." << std::endl;
283 aggStyle =
"Point Cloud";
284 this->doPointCloud(vertices, geomSize, numLocalAggs, numFineNodes);
289 if (pL.
get<
bool>(
"visualization: fine graph edges")) {
291 "Could not get information about fine graph.");
293 std::vector<LocalOrdinal> fine_edges_vertices;
294 std::vector<LocalOrdinal> fine_edges_geomSize;
295 this->doGraphEdges(fine_edges_vertices, fine_edges_geomSize, fineGraph, xCoords, yCoords, zCoords);
297 std::string fEdgeFineFile = this->getFileName(comm->getSize(), comm->getRank(), fineLevel.GetLevelID(), pL);
298 std::string fEdgeFile = fEdgeFineFile.insert(fEdgeFineFile.rfind(
".vtu"),
"-finegraph");
299 std::ofstream edgeStream(fEdgeFile.c_str());
301 std::vector<int> uniqueFineEdges = this->makeUnique(fine_edges_vertices);
302 this->writeFileVTKOpening(edgeStream, uniqueFineEdges, fine_edges_geomSize);
303 this->writeFileVTKNodes(edgeStream, uniqueFineEdges, nodeMap);
304 this->writeFileVTKData(edgeStream, uniqueFineEdges, myAggOffset, vertex2AggId, comm->getRank());
305 this->writeFileVTKCoordinates(edgeStream, uniqueFineEdges, xCoords, yCoords, zCoords, coords->getNumVectors());
306 this->writeFileVTKCells(edgeStream, uniqueFineEdges, fine_edges_vertices, fine_edges_geomSize);
307 this->writeFileVTKClosing(edgeStream);
312 #if 0 // we don't have access to the coarse graph
313 if (pL.
get<
bool>(
"visualization: coarse graph edges")) {
314 RCP<LWGraph> coarseGraph = Get<RCP<LWGraph> >(coarseLevel,
"Graph");
320 coarseghostedCoords->doImport(*coarsecoords, *coarsecoordImporter,
Xpetra::INSERT);
321 coarsecoords = coarseghostedCoords;
326 if(coarsecoords->getNumVectors() == 3) {
332 std::vector<LocalOrdinal> coarse_edges_vertices;
333 std::vector<LocalOrdinal> coarse_edges_geomSize;
334 this->doGraphEdges(coarse_edges_vertices, coarse_edges_geomSize, coarseGraph, cx, cy, cz);
336 std::string cEdgeFineFile = this->getFileName(comm->getSize(), comm->getRank(), coarseLevel.GetLevelID(), pL);
337 std::string cEdgeFile = cEdgeFineFile.insert(cEdgeFineFile.rfind(
".vtu"),
"-coarsegraph");
338 std::ofstream cedgeStream(cEdgeFile.c_str());
340 std::vector<int> uniqueCoarseEdges = this->makeUnique(coarse_edges_vertices);
341 this->writeFileVTKOpening(cedgeStream, uniqueCoarseEdges, coarse_edges_geomSize);
342 this->writeFileVTKNodes(cedgeStream, uniqueCoarseEdges, coarsenodeMap);
344 this->writeFileVTKCoordinates(cedgeStream, uniqueCoarseEdges, cx, cy, cz, coarsecoords->getNumVectors());
345 this->writeFileVTKCells(cedgeStream, uniqueCoarseEdges, coarse_edges_vertices, coarse_edges_geomSize);
346 this->writeFileVTKClosing(cedgeStream);
351 if (pL.
get<
int>(
"visualization: start level") <= fineLevel.GetLevelID()) {
353 std::string filenameToWrite = this->getFileName(comm->getSize(), comm->getRank(), fineLevel.GetLevelID(), pL);
354 std::ofstream fout(filenameToWrite.c_str());
356 std::vector<int> uniqueFine = this->makeUnique(vertices);
357 this->writeFileVTKOpening(fout, uniqueFine, geomSize);
358 this->writeFileVTKNodes(fout, uniqueFine, nodeMap);
359 this->writeFileVTKData(fout, uniqueFine, myAggOffset, vertex2AggId, comm->getRank());
360 this->writeFileVTKCoordinates(fout, uniqueFine, xCoords, yCoords, zCoords, coords->getNumVectors());
361 this->writeFileVTKCells(fout, uniqueFine, vertices, geomSize);
362 this->writeFileVTKClosing(fout);
366 if (comm->getRank() == 0) {
367 std::string pvtuFilename = this->getPVTUFileName(comm->getSize(), comm->getRank(), fineLevel.GetLevelID(), pL);
368 std::string baseFname = this->getBaseFileName(comm->getSize(), fineLevel.GetLevelID(), pL);
369 std::ofstream pvtu(pvtuFilename.c_str());
370 this->writePVTU(pvtu, baseFname, comm->getSize(), pL.
get<
bool>(
"visualization: fine graph edges"));
375 if (comm->getRank() == 0 && pL.
get<
bool>(
"visualization: build colormap")) {
376 this->buildColormap();
Important warning messages (one line)
MueLu::DefaultLocalOrdinal LocalOrdinal
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)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
RCP< ParameterList > GetValidParameterList() const
bool isParameter(const std::string &name) const
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Class that holds all level-specific information.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
const RCP< const Map > GetImportMap() const
Return overlapping import map (nodes).
Exception throws to report errors in the internal logical of the program.
Description of what is happening (more verbose)
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
static RCP< Import< LocalOrdinal, GlobalOrdinal, Node > > Build(const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &source, const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &target, const Teuchos::RCP< Teuchos::ParameterList > &plist=Teuchos::null)