47#ifndef MUELU_COARSENINGVISUALIZATIONFACTORY_DEF_HPP_
48#define MUELU_COARSENINGVISUALIZATIONFACTORY_DEF_HPP_
50#include <Xpetra_Matrix.hpp>
51#include <Xpetra_ImportFactory.hpp>
52#include <Xpetra_MultiVectorFactory.hpp>
59 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
64 validParamList->set<
int > (
"visualization: start level", 0,
"visualize only levels with level ids greater or equal than start level");
66 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.");
67 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");
68 validParamList->set< RCP<const FactoryBase> >(
"Coordinates", Teuchos::null,
"Factory for Coordinates.");
69 validParamList->set< RCP<const FactoryBase> >(
"Graph", Teuchos::null,
"Factory for Graph.");
71 return validParamList;
74 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
76 this->Input(fineLevel,
"Coordinates");
78 const ParameterList & pL = this->GetParameterList();
79 TEUCHOS_TEST_FOR_EXCEPTION(pL.isParameter(
"P") && GetFactory(
"P") != Teuchos::null &&
81 "You must not declare both P and Ptent. Use only once for visualization.");
82 TEUCHOS_TEST_FOR_EXCEPTION(GetFactory(
"P") == Teuchos::null && GetFactory(
"Ptent") == Teuchos::null,
Exceptions::RuntimeError,
83 "You have to either declare P or Ptent for visualization, but not both.");
85 if (GetFactory(
"P") != Teuchos::null && GetFactory(
"Ptent") == Teuchos::null)
86 this->Input(coarseLevel,
"P");
87 else if (GetFactory(
"Ptent") != Teuchos::null && GetFactory(
"P") == Teuchos::null)
88 this->Input(coarseLevel,
"Ptent");
90 if(pL.get<
bool>(
"visualization: fine graph edges"))
91 Input(fineLevel,
"Graph");
93 if(pL.get<
bool>(
"visualization: coarse graph edges")) {
94 Input(coarseLevel,
"Coordinates");
95 Input(coarseLevel,
"Graph");
100 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
103 RCP<GraphBase> fineGraph = Teuchos::null;
104 RCP<Matrix> P = Teuchos::null;
105 const ParameterList & pL = this->GetParameterList();
106 if (this->GetFactory(
"P") != Teuchos::null && this->GetFactory(
"Ptent") == Teuchos::null)
107 P = Get< RCP<Matrix> >(coarseLevel,
"P");
108 if (GetFactory(
"Ptent") != Teuchos::null && GetFactory(
"P") == Teuchos::null)
109 P = Get< RCP<Matrix> >(coarseLevel,
"Ptent");
111 RCP<const Teuchos::Comm<int> > comm = P->getRowMap()->getComm();
115 RCP<const StridedMap> strRowMap = Teuchos::null;
116 if (P->IsView(
"stridedMaps") && Teuchos::rcp_dynamic_cast<const StridedMap>(P->getRowMap(
"stridedMaps")) != Teuchos::null) {
117 strRowMap = Teuchos::rcp_dynamic_cast<const StridedMap>(P->getRowMap(
"stridedMaps"));
120 std::vector<size_t> stridingInfo = strRowMap->getStridingData();
121 for (
size_t j = 0; j < Teuchos::as<size_t>(blockid); j++)
122 stridedRowOffset += stridingInfo[j];
123 dofsPerNode = Teuchos::as<LocalOrdinal>(stridingInfo[blockid]);
125 dofsPerNode = strRowMap->getFixedBlockSize();
127 GetOStream(
Runtime1) <<
"CoarseningVisualizationFactory::Build():" <<
" #dofs per node = " << dofsPerNode << std::endl;
132 RCP<const StridedMap> strDomainMap = Teuchos::null;
133 if (P->IsView(
"stridedMaps") && Teuchos::rcp_dynamic_cast<const StridedMap>(P->getRowMap(
"stridedMaps")) != Teuchos::null) {
134 strDomainMap = Teuchos::rcp_dynamic_cast<const StridedMap>(P->getColMap(
"stridedMaps"));
135 LocalOrdinal blockid = strDomainMap->getStridedBlockId();
138 std::vector<size_t> stridingInfo = strDomainMap->getStridingData();
139 for (
size_t j = 0; j < Teuchos::as<size_t>(blockid); j++)
140 stridedColumnOffset += stridingInfo[j];
141 columnsPerNode = Teuchos::as<LocalOrdinal>(stridingInfo[blockid]);
143 columnsPerNode = strDomainMap->getFixedBlockSize();
145 GetOStream(
Runtime1) <<
"CoarseningVisualizationFactory::Build():" <<
" #columns per node = " << columnsPerNode << std::endl;
149 TEUCHOS_TEST_FOR_EXCEPTION(strDomainMap->getLocalNumElements() != P->getColMap()->getLocalNumElements(),
Exceptions::RuntimeError,
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) {
158 ArrayView<const LO> indices;
159 ArrayView<const SC> vals;
160 P->getLocalRowView(row, indices, vals);
162 for(
typename ArrayView<const LO>::iterator c = indices.begin(); c != indices.end(); ++c) {
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];
190 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::as<LO>(P->getRowMap()->getLocalNumElements()) / dofsPerNode != Teuchos::as<LocalOrdinal>(coords->getLocalLength()),
Exceptions::RuntimeError,
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<GraphBase> >(fineLevel,
"Graph");
197 RCP<Import> coordImporter = Xpetra::ImportFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(coords->getMap(), fineGraph->GetImportMap());
198 RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType,
LocalOrdinal,
GlobalOrdinal,
Node> > ghostedCoords = Xpetra::MultiVectorFactory<typename Teuchos::ScalarTraits<Scalar>::coordinateType,
LocalOrdinal,
GlobalOrdinal,
Node>::Build(fineGraph->GetImportMap(), coords->getNumVectors());
199 ghostedCoords->doImport(*coords, *coordImporter, Xpetra::INSERT);
200 coords = ghostedCoords;
203 Teuchos::RCP<const Map> nodeMap = coords->getMap();
205 Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> xCoords = Teuchos::arcp_reinterpret_cast<const typename Teuchos::ScalarTraits<Scalar>::coordinateType>(coords->getData(0));
206 Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> yCoords = Teuchos::arcp_reinterpret_cast<const typename Teuchos::ScalarTraits<Scalar>::coordinateType>(coords->getData(1));
207 Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> zCoords = Teuchos::null;
208 if(coords->getNumVectors() == 3) {
209 zCoords = Teuchos::arcp_reinterpret_cast<const typename Teuchos::ScalarTraits<Scalar>::coordinateType>(coords->getData(2));
213 LocalOrdinal numFineNodes = Teuchos::as<LocalOrdinal>(coords->getLocalLength());
216 ArrayRCP<LocalOrdinal> vertex2AggId(numFineNodes, -1);
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);
230 typename Teuchos::ScalarTraits<Scalar>::coordinateType xCenter = 0.0;
231 typename Teuchos::ScalarTraits<Scalar>::coordinateType yCenter = 0.0;
232 typename Teuchos::ScalarTraits<Scalar>::coordinateType zCenter = 0.0;
235 for(
typename std::set<LocalOrdinal>::iterator it = localAggs[i].begin(); it != localAggs[i].end(); ++it) {
236 xCenter += xCoords[*it];
237 yCenter += yCoords[*it];
238 if(coords->getNumVectors() == 3) zCenter += zCoords[*it];
240 xCenter /= Teuchos::as<LocalOrdinal>(localAggs[i].size());
241 yCenter /= Teuchos::as<LocalOrdinal>(localAggs[i].size());
242 zCenter /= Teuchos::as<LocalOrdinal>(localAggs[i].size());
246 typename Teuchos::ScalarTraits<Scalar>::coordinateType minDistance = Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<Scalar>::coordinateType>::one() / Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<Scalar>::coordinateType>::sfmin();
247 for(
typename std::set<LocalOrdinal>::iterator it = localAggs[i].begin(); it != localAggs[i].end(); ++it) {
248 typename Teuchos::ScalarTraits<Scalar>::coordinateType tempx = xCenter - xCoords[*it];
249 typename Teuchos::ScalarTraits<Scalar>::coordinateType tempy = yCenter - yCoords[*it];
250 typename Teuchos::ScalarTraits<Scalar>::coordinateType tempz = 0.0;
251 if(coords->getNumVectors() == 3) tempz = zCenter - zCoords[*it];
252 typename Teuchos::ScalarTraits<Scalar>::coordinateType mydistance = 0.0;
253 mydistance += tempx*tempx;
254 mydistance += tempy*tempy;
255 mydistance += tempz*tempz;
256 mydistance = sqrt(mydistance);
257 if (mydistance <= minDistance) {
258 minDistance = mydistance;
263 isRoot[rootCandidate] =
true;
266 std::vector<LocalOrdinal> vertices;
267 std::vector<LocalOrdinal> geomSize;
268 int vizLevel = pL.get<
int>(
"visualization: start level");
271 std::string aggStyle = pL.get<std::string>(
"visualization: style");
272 if(aggStyle ==
"Point Cloud")
273 this->doPointCloud(vertices, geomSize, numLocalAggs, numFineNodes);
274 else if(aggStyle ==
"Jacks")
275 this->doJacks(vertices, geomSize, numLocalAggs, numFineNodes, isRoot, vertex2AggId);
276 else if(aggStyle ==
"Convex Hulls") {
279 if(coords->getNumVectors() == 3)
280 this->doConvexHulls3D(vertices, geomSize, numLocalAggs, numFineNodes, isRoot, vertex2AggId, xCoords, yCoords, zCoords);
281 else if(coords->getNumVectors() == 2)
282 this->doConvexHulls2D(vertices, geomSize, numLocalAggs, numFineNodes, isRoot, vertex2AggId, xCoords, yCoords);
284 else if(aggStyle ==
"CGAL Convex Hulls") {
285#ifdef HAVE_MUELU_CGAL
286 if(coords->getNumVectors() == 3)
287 this->doCGALConvexHulls3D(vertices, geomSize, numLocalAggs, numFineNodes, isRoot, vertex2AggId, xCoords, yCoords, zCoords);
288 else if(coords->getNumVectors() == 2)
289 this->doCGALConvexHulls2D(vertices, geomSize, numLocalAggs, numFineNodes, isRoot, vertex2AggId, xCoords, yCoords);
294 GetOStream(
Warnings0) <<
" Warning: Unrecognized agg style.\nPossible values are Point Cloud, Jacks, Convex Hulls.\nDefaulting to Point Cloud." << std::endl;
295 aggStyle =
"Point Cloud";
296 this->doPointCloud(vertices, geomSize, numLocalAggs, numFineNodes);
301 if(pL.get<
bool>(
"visualization: fine graph edges")) {
303 "Could not get information about fine graph.");
305 std::vector<LocalOrdinal> fine_edges_vertices;
306 std::vector<LocalOrdinal> fine_edges_geomSize;
307 this->doGraphEdges(fine_edges_vertices, fine_edges_geomSize, fineGraph, xCoords, yCoords, zCoords);
309 std::string fEdgeFineFile = this->getFileName(comm->getSize(), comm->getRank(), fineLevel.
GetLevelID(), pL);
310 std::string fEdgeFile = fEdgeFineFile.insert(fEdgeFineFile.rfind(
".vtu"),
"-finegraph");
311 std::ofstream edgeStream(fEdgeFile.c_str());
313 std::vector<int> uniqueFineEdges = this->makeUnique(fine_edges_vertices);
314 this->writeFileVTKOpening(edgeStream, uniqueFineEdges, fine_edges_geomSize);
315 this->writeFileVTKNodes(edgeStream, uniqueFineEdges, nodeMap);
316 this->writeFileVTKData(edgeStream, uniqueFineEdges, myAggOffset, vertex2AggId, comm->getRank());
317 this->writeFileVTKCoordinates(edgeStream, uniqueFineEdges, xCoords, yCoords, zCoords, coords->getNumVectors());
318 this->writeFileVTKCells(edgeStream, uniqueFineEdges, fine_edges_vertices, fine_edges_geomSize);
319 this->writeFileVTKClosing(edgeStream);
325 if (pL.get<
bool>(
"visualization: coarse graph edges")) {
326 RCP<GraphBase> coarseGraph = Get<RCP<GraphBase> >(coarseLevel,
"Graph");
328 Teuchos::RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType,
LocalOrdinal,
GlobalOrdinal,
Node> > coarsecoords = Get<RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType,
LocalOrdinal,
GlobalOrdinal,
Node> > >(coarseLevel,
"Coordinates");
330 RCP<Import> coarsecoordImporter = Xpetra::ImportFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(coarsecoords->getMap(), coarseGraph->GetImportMap());
331 RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType,
LocalOrdinal,
GlobalOrdinal,
Node> > coarseghostedCoords = Xpetra::MultiVectorFactory<typename Teuchos::ScalarTraits<Scalar>::coordinateType,
LocalOrdinal,
GlobalOrdinal,
Node>::Build(coarseGraph->GetImportMap(), coarsecoords->getNumVectors());
332 coarseghostedCoords->doImport(*coarsecoords, *coarsecoordImporter, Xpetra::INSERT);
333 coarsecoords = coarseghostedCoords;
335 Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> cx = Teuchos::arcp_reinterpret_cast<const typename Teuchos::ScalarTraits<Scalar>::coordinateType>(coarsecoords->getData(0));
336 Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> cy = Teuchos::arcp_reinterpret_cast<const typename Teuchos::ScalarTraits<Scalar>::coordinateType>(coarsecoords->getData(1));
337 Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> cz = Teuchos::null;
338 if(coarsecoords->getNumVectors() == 3) {
339 cz = Teuchos::arcp_reinterpret_cast<const typename Teuchos::ScalarTraits<Scalar>::coordinateType>(coarsecoords->getData(2));
342 Teuchos::RCP<const Map> coarsenodeMap = coarsecoords->getMap();
344 std::vector<LocalOrdinal> coarse_edges_vertices;
345 std::vector<LocalOrdinal> coarse_edges_geomSize;
346 this->doGraphEdges(coarse_edges_vertices, coarse_edges_geomSize, coarseGraph, cx, cy, cz);
348 std::string cEdgeFineFile = this->getFileName(comm->getSize(), comm->getRank(), coarseLevel.GetLevelID(), pL);
349 std::string cEdgeFile = cEdgeFineFile.insert(cEdgeFineFile.rfind(
".vtu"),
"-coarsegraph");
350 std::ofstream cedgeStream(cEdgeFile.c_str());
352 std::vector<int> uniqueCoarseEdges = this->makeUnique(coarse_edges_vertices);
353 this->writeFileVTKOpening(cedgeStream, uniqueCoarseEdges, coarse_edges_geomSize);
354 this->writeFileVTKNodes(cedgeStream, uniqueCoarseEdges, coarsenodeMap);
356 this->writeFileVTKCoordinates(cedgeStream, uniqueCoarseEdges, cx, cy, cz, coarsecoords->getNumVectors());
357 this->writeFileVTKCells(cedgeStream, uniqueCoarseEdges, coarse_edges_vertices, coarse_edges_geomSize);
358 this->writeFileVTKClosing(cedgeStream);
363 if(pL.get<
int>(
"visualization: start level") <= fineLevel.
GetLevelID()) {
365 std::string filenameToWrite = this->getFileName(comm->getSize(), comm->getRank(), fineLevel.
GetLevelID(), pL);
366 std::ofstream fout (filenameToWrite.c_str());
368 std::vector<int> uniqueFine = this->makeUnique(vertices);
369 this->writeFileVTKOpening(fout, uniqueFine, geomSize);
370 this->writeFileVTKNodes(fout, uniqueFine, nodeMap);
371 this->writeFileVTKData(fout, uniqueFine, myAggOffset, vertex2AggId, comm->getRank());
372 this->writeFileVTKCoordinates(fout, uniqueFine, xCoords, yCoords, zCoords, coords->getNumVectors());
373 this->writeFileVTKCells(fout, uniqueFine, vertices, geomSize);
374 this->writeFileVTKClosing(fout);
378 if(comm->getRank() == 0) {
379 std::string pvtuFilename = this->getPVTUFileName(comm->getSize(), comm->getRank(), fineLevel.
GetLevelID(), pL);
380 std::string baseFname = this->getBaseFileName(comm->getSize(), fineLevel.
GetLevelID(), pL);
381 std::ofstream pvtu(pvtuFilename.c_str());
382 this->writePVTU(pvtu, baseFname, comm->getSize(), pL.get<
bool>(
"visualization: fine graph edges"));
387 if(comm->getRank() == 0 && pL.get<
bool>(
"visualization: build colormap")) {
388 this->buildColormap();
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultGlobalOrdinal GlobalOrdinal
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
Exception throws to report errors in the internal logical of the program.
Class that holds all level-specific information.
int GetLevelID() const
Return level number.
RCP< ParameterList > GetValidParameterList() const
Namespace for MueLu classes and methods.
@ Warnings0
Important warning messages (one line)
@ Runtime1
Description of what is happening (more verbose)