MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_CoarseningVisualizationFactory_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // MueLu: A package for multigrid based preconditioning
4 //
5 // Copyright 2012 NTESS and the MueLu contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef MUELU_COARSENINGVISUALIZATIONFACTORY_DEF_HPP_
11 #define MUELU_COARSENINGVISUALIZATIONFACTORY_DEF_HPP_
12 
13 #include <Xpetra_Matrix.hpp>
14 #include <Xpetra_ImportFactory.hpp>
15 #include <Xpetra_MultiVectorFactory.hpp>
17 #include "MueLu_Level.hpp"
18 
19 namespace MueLu {
20 
21 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
24 
25  validParamList->set<int>("visualization: start level", 0, "visualize only levels with level ids greater or equal than start level"); // Remove me?
26 
27  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.");
28  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");
29  validParamList->set<RCP<const FactoryBase> >("Coordinates", Teuchos::null, "Factory for Coordinates.");
30  validParamList->set<RCP<const FactoryBase> >("Graph", Teuchos::null, "Factory for Graph.");
31 
32  return validParamList;
33 }
34 
35 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
37  this->Input(fineLevel, "Coordinates");
38 
39  const ParameterList &pL = this->GetParameterList();
40  TEUCHOS_TEST_FOR_EXCEPTION(pL.isParameter("P") && GetFactory("P") != Teuchos::null &&
41  pL.isParameter("Ptent") && GetFactory("Ptent") != Teuchos::null,
43  "You must not declare both P and Ptent. Use only once for visualization.");
44  TEUCHOS_TEST_FOR_EXCEPTION(GetFactory("P") == Teuchos::null && GetFactory("Ptent") == Teuchos::null, Exceptions::RuntimeError,
45  "You have to either declare P or Ptent for visualization, but not both.");
46 
47  if (GetFactory("P") != Teuchos::null && GetFactory("Ptent") == Teuchos::null)
48  this->Input(coarseLevel, "P");
49  else if (GetFactory("Ptent") != Teuchos::null && GetFactory("P") == Teuchos::null)
50  this->Input(coarseLevel, "Ptent");
51 
52  if (pL.get<bool>("visualization: fine graph edges"))
53  Input(fineLevel, "Graph");
54 #if 0
55  if(pL.get<bool>("visualization: coarse graph edges")) {
56  Input(coarseLevel, "Coordinates");
57  Input(coarseLevel, "Graph");
58  }
59 #endif
60 }
61 
62 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
64  RCP<LWGraph> fineGraph = Teuchos::null;
65  RCP<Matrix> P = Teuchos::null;
66  const ParameterList &pL = this->GetParameterList();
67  if (this->GetFactory("P") != Teuchos::null && this->GetFactory("Ptent") == Teuchos::null)
68  P = Get<RCP<Matrix> >(coarseLevel, "P");
69  if (GetFactory("Ptent") != Teuchos::null && GetFactory("P") == Teuchos::null)
70  P = Get<RCP<Matrix> >(coarseLevel, "Ptent");
71 
72  RCP<const Teuchos::Comm<int> > comm = P->getRowMap()->getComm();
73 
74  LocalOrdinal dofsPerNode = 1;
75  LocalOrdinal stridedRowOffset = 0;
76  RCP<const StridedMap> strRowMap = Teuchos::null;
77  if (P->IsView("stridedMaps") && Teuchos::rcp_dynamic_cast<const StridedMap>(P->getRowMap("stridedMaps")) != Teuchos::null) {
78  strRowMap = Teuchos::rcp_dynamic_cast<const StridedMap>(P->getRowMap("stridedMaps"));
79  LocalOrdinal blockid = strRowMap->getStridedBlockId();
80  if (blockid > -1) {
81  std::vector<size_t> stridingInfo = strRowMap->getStridingData();
82  for (size_t j = 0; j < Teuchos::as<size_t>(blockid); j++)
83  stridedRowOffset += stridingInfo[j];
84  dofsPerNode = Teuchos::as<LocalOrdinal>(stridingInfo[blockid]);
85  } else {
86  dofsPerNode = strRowMap->getFixedBlockSize();
87  }
88  GetOStream(Runtime1) << "CoarseningVisualizationFactory::Build():"
89  << " #dofs per node = " << dofsPerNode << std::endl;
90  }
91 
92  LocalOrdinal columnsPerNode = dofsPerNode;
93  LocalOrdinal stridedColumnOffset = 0;
94  RCP<const StridedMap> strDomainMap = Teuchos::null;
95  if (P->IsView("stridedMaps") && Teuchos::rcp_dynamic_cast<const StridedMap>(P->getRowMap("stridedMaps")) != Teuchos::null) {
96  strDomainMap = Teuchos::rcp_dynamic_cast<const StridedMap>(P->getColMap("stridedMaps"));
97  LocalOrdinal blockid = strDomainMap->getStridedBlockId();
98 
99  if (blockid > -1) {
100  std::vector<size_t> stridingInfo = strDomainMap->getStridingData();
101  for (size_t j = 0; j < Teuchos::as<size_t>(blockid); j++)
102  stridedColumnOffset += stridingInfo[j];
103  columnsPerNode = Teuchos::as<LocalOrdinal>(stridingInfo[blockid]);
104  } else {
105  columnsPerNode = strDomainMap->getFixedBlockSize();
106  }
107  GetOStream(Runtime1) << "CoarseningVisualizationFactory::Build():"
108  << " #columns per node = " << columnsPerNode << std::endl;
109  }
110 
111  // TODO add support for overlapping aggregates
112  TEUCHOS_TEST_FOR_EXCEPTION(strDomainMap->getLocalNumElements() != P->getColMap()->getLocalNumElements(), Exceptions::RuntimeError,
113  "CoarseningVisualization only supports non-overlapping transfers");
114 
115  // number of local "aggregates"
116  LocalOrdinal numLocalAggs = strDomainMap->getLocalNumElements() / columnsPerNode;
117  std::vector<std::set<LocalOrdinal> > localAggs(numLocalAggs);
118 
119  // do loop over all local rows of prolongator and extract column information
120  for (LO row = 0; row < Teuchos::as<LO>(P->getRowMap()->getLocalNumElements()); ++row) {
121  ArrayView<const LO> indices;
122  ArrayView<const SC> vals;
123  P->getLocalRowView(row, indices, vals);
124 
125  for (typename ArrayView<const LO>::iterator c = indices.begin(); c != indices.end(); ++c) {
126  localAggs[(*c) / columnsPerNode].insert(row / dofsPerNode);
127  }
128  }
129 
130  // determine number of "aggs" per proc and calculate local "agg" offset...
131  std::vector<int> myLocalAggsPerProc(comm->getSize(), 0);
132  std::vector<int> numLocalAggsPerProc(comm->getSize(), 0);
133  myLocalAggsPerProc[comm->getRank()] = numLocalAggs;
134  Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, comm->getSize(), &myLocalAggsPerProc[0], &numLocalAggsPerProc[0]);
135 
136  LocalOrdinal myAggOffset = 0;
137  for (int i = 0; i < comm->getRank(); ++i) {
138  myAggOffset += numLocalAggsPerProc[i];
139  }
140 
141  /*for (LocalOrdinal i = 0; i < numLocalAggs; ++i) {
142 
143  std::cout << "PROC: " << comm->getRank() << " Local aggregate: " << i + myAggOffset << " with nodes: ";
144  for( typename std::set<LocalOrdinal>::iterator it = localAggs[i].begin(); it != localAggs[i].end(); ++it) {
145  std::cout << *it << ", ";
146  }
147  std::cout << std::endl;
148  }*/
149 
150  // get fine level coordinate information
151  Teuchos::RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType, LocalOrdinal, GlobalOrdinal, Node> > coords = Get<RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType, LocalOrdinal, GlobalOrdinal, Node> > >(fineLevel, "Coordinates");
152 
153  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::as<LO>(P->getRowMap()->getLocalNumElements()) / dofsPerNode != Teuchos::as<LocalOrdinal>(coords->getLocalLength()), Exceptions::RuntimeError,
154  "Number of fine level nodes in coordinates is inconsistent with dof based information");
155 
156  // communicate fine level coordinates if necessary
157  if (pL.get<bool>("visualization: fine graph edges")) {
158  fineGraph = Get<RCP<LWGraph> >(fineLevel, "Graph");
159 
160  RCP<Import> coordImporter = Xpetra::ImportFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(coords->getMap(), fineGraph->GetImportMap());
162  ghostedCoords->doImport(*coords, *coordImporter, Xpetra::INSERT);
163  coords = ghostedCoords;
164  }
165 
166  Teuchos::RCP<const Map> nodeMap = coords->getMap();
167 
168  Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> xCoords = Teuchos::arcp_reinterpret_cast<const typename Teuchos::ScalarTraits<Scalar>::coordinateType>(coords->getData(0));
169  Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> yCoords = Teuchos::arcp_reinterpret_cast<const typename Teuchos::ScalarTraits<Scalar>::coordinateType>(coords->getData(1));
171  if (coords->getNumVectors() == 3) {
172  zCoords = Teuchos::arcp_reinterpret_cast<const typename Teuchos::ScalarTraits<Scalar>::coordinateType>(coords->getData(2));
173  }
174 
175  // determine number of nodes on fine level
176  LocalOrdinal numFineNodes = Teuchos::as<LocalOrdinal>(coords->getLocalLength());
177 
178  // create vertex2AggId array
179  ArrayRCP<LocalOrdinal> vertex2AggId(numFineNodes, -1);
180  for (LocalOrdinal i = 0; i < numLocalAggs; ++i) {
181  // TODO: check if entry = -1
182  for (typename std::set<LocalOrdinal>::iterator it = localAggs[i].begin(); it != localAggs[i].end(); ++it) {
183  vertex2AggId[*it] = i;
184  }
185  }
186 
187  // we have no information which node is root and which is not
188  // we could either look at the entries in P again or build some new root nodes
189  // assuming that root nodes are usually at the center of the aggregate
190  std::vector<bool> isRoot(numFineNodes, false);
191  for (LocalOrdinal i = 0; i < numLocalAggs; ++i) {
192  typename Teuchos::ScalarTraits<Scalar>::coordinateType xCenter = 0.0;
193  typename Teuchos::ScalarTraits<Scalar>::coordinateType yCenter = 0.0;
194  typename Teuchos::ScalarTraits<Scalar>::coordinateType zCenter = 0.0;
195 
196  // loop over all nodes in aggregate i and determine center coordinates of aggregate
197  for (typename std::set<LocalOrdinal>::iterator it = localAggs[i].begin(); it != localAggs[i].end(); ++it) {
198  xCenter += xCoords[*it];
199  yCenter += yCoords[*it];
200  if (coords->getNumVectors() == 3) zCenter += zCoords[*it];
201  }
202  xCenter /= Teuchos::as<LocalOrdinal>(localAggs[i].size());
203  yCenter /= Teuchos::as<LocalOrdinal>(localAggs[i].size());
204  zCenter /= Teuchos::as<LocalOrdinal>(localAggs[i].size());
205 
206  // loop over all nodes in aggregate i and find node which is closest to aggregate center
207  LocalOrdinal rootCandidate = -1;
209  for (typename std::set<LocalOrdinal>::iterator it = localAggs[i].begin(); it != localAggs[i].end(); ++it) {
210  typename Teuchos::ScalarTraits<Scalar>::coordinateType tempx = xCenter - xCoords[*it];
211  typename Teuchos::ScalarTraits<Scalar>::coordinateType tempy = yCenter - yCoords[*it];
213  if (coords->getNumVectors() == 3) tempz = zCenter - zCoords[*it];
214  typename Teuchos::ScalarTraits<Scalar>::coordinateType mydistance = 0.0;
215  mydistance += tempx * tempx;
216  mydistance += tempy * tempy;
217  mydistance += tempz * tempz;
218  mydistance = sqrt(mydistance);
219  if (mydistance <= minDistance) {
220  minDistance = mydistance;
221  rootCandidate = *it;
222  }
223  }
224 
225  isRoot[rootCandidate] = true;
226  }
227 
228  std::vector<LocalOrdinal> vertices;
229  std::vector<LocalOrdinal> geomSize;
230  int vizLevel = pL.get<int>("visualization: start level");
231  if (vizLevel <= fineLevel.GetLevelID()) {
232  std::string aggStyle = pL.get<std::string>("visualization: style");
233  if (aggStyle == "Point Cloud")
234  this->doPointCloud(vertices, geomSize, numLocalAggs, numFineNodes);
235  else if (aggStyle == "Jacks")
236  this->doJacks(vertices, geomSize, numLocalAggs, numFineNodes, isRoot, vertex2AggId);
237  else if (aggStyle == "Convex Hulls") {
238  // TODO do a smarter distinction and check the number of z levels...
239  // loop over all coordinates and collect coordinate components in sets...
240  if (coords->getNumVectors() == 3)
241  this->doConvexHulls3D(vertices, geomSize, numLocalAggs, numFineNodes, isRoot, vertex2AggId, xCoords, yCoords, zCoords);
242  else if (coords->getNumVectors() == 2)
243  this->doConvexHulls2D(vertices, geomSize, numLocalAggs, numFineNodes, isRoot, vertex2AggId, xCoords, yCoords);
244  } else {
245  GetOStream(Warnings0) << " Warning: Unrecognized agg style.\nPossible values are Point Cloud, Jacks, Convex Hulls.\nDefaulting to Point Cloud." << std::endl;
246  aggStyle = "Point Cloud";
247  this->doPointCloud(vertices, geomSize, numLocalAggs, numFineNodes);
248  }
249  }
250 
251  // write out fine edge information
252  if (pL.get<bool>("visualization: fine graph edges")) {
253  TEUCHOS_TEST_FOR_EXCEPTION(fineGraph == Teuchos::null, Exceptions::RuntimeError,
254  "Could not get information about fine graph.");
255 
256  std::vector<LocalOrdinal> fine_edges_vertices;
257  std::vector<LocalOrdinal> fine_edges_geomSize;
258  this->doGraphEdges(fine_edges_vertices, fine_edges_geomSize, fineGraph, xCoords, yCoords, zCoords);
259 
260  std::string fEdgeFineFile = this->getFileName(comm->getSize(), comm->getRank(), fineLevel.GetLevelID(), pL);
261  std::string fEdgeFile = fEdgeFineFile.insert(fEdgeFineFile.rfind(".vtu"), "-finegraph");
262  std::ofstream edgeStream(fEdgeFile.c_str());
263 
264  std::vector<int> uniqueFineEdges = this->makeUnique(fine_edges_vertices);
265  this->writeFileVTKOpening(edgeStream, uniqueFineEdges, fine_edges_geomSize);
266  this->writeFileVTKNodes(edgeStream, uniqueFineEdges, nodeMap);
267  this->writeFileVTKData(edgeStream, uniqueFineEdges, myAggOffset, vertex2AggId, comm->getRank());
268  this->writeFileVTKCoordinates(edgeStream, uniqueFineEdges, xCoords, yCoords, zCoords, coords->getNumVectors());
269  this->writeFileVTKCells(edgeStream, uniqueFineEdges, fine_edges_vertices, fine_edges_geomSize);
270  this->writeFileVTKClosing(edgeStream);
271  edgeStream.close();
272  }
273 
274  // communicate fine level coordinates if necessary
275 #if 0 // we don't have access to the coarse graph
276  if (pL.get<bool>("visualization: coarse graph edges")) {
277  RCP<LWGraph> coarseGraph = Get<RCP<LWGraph> >(coarseLevel, "Graph");
278 
279  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");
280 
281  RCP<Import> coarsecoordImporter = Xpetra::ImportFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(coarsecoords->getMap(), coarseGraph->GetImportMap());
283  coarseghostedCoords->doImport(*coarsecoords, *coarsecoordImporter, Xpetra::INSERT);
284  coarsecoords = coarseghostedCoords;
285 
286  Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> cx = Teuchos::arcp_reinterpret_cast<const typename Teuchos::ScalarTraits<Scalar>::coordinateType>(coarsecoords->getData(0));
287  Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> cy = Teuchos::arcp_reinterpret_cast<const typename Teuchos::ScalarTraits<Scalar>::coordinateType>(coarsecoords->getData(1));
289  if(coarsecoords->getNumVectors() == 3) {
290  cz = Teuchos::arcp_reinterpret_cast<const typename Teuchos::ScalarTraits<Scalar>::coordinateType>(coarsecoords->getData(2));
291  }
292 
293  Teuchos::RCP<const Map> coarsenodeMap = coarsecoords->getMap();
294 
295  std::vector<LocalOrdinal> coarse_edges_vertices;
296  std::vector<LocalOrdinal> coarse_edges_geomSize;
297  this->doGraphEdges(coarse_edges_vertices, coarse_edges_geomSize, coarseGraph, cx, cy, cz);
298 
299  std::string cEdgeFineFile = this->getFileName(comm->getSize(), comm->getRank(), coarseLevel.GetLevelID(), pL);
300  std::string cEdgeFile = cEdgeFineFile.insert(cEdgeFineFile.rfind(".vtu"), "-coarsegraph");
301  std::ofstream cedgeStream(cEdgeFile.c_str());
302 
303  std::vector<int> uniqueCoarseEdges = this->makeUnique(coarse_edges_vertices);
304  this->writeFileVTKOpening(cedgeStream, uniqueCoarseEdges, coarse_edges_geomSize);
305  this->writeFileVTKNodes(cedgeStream, uniqueCoarseEdges, coarsenodeMap);
306  //this->writeFileVTKData(edgeStream, uniqueCoarseEdges, myAggOffset, vertex2AggId, comm->getRank());
307  this->writeFileVTKCoordinates(cedgeStream, uniqueCoarseEdges, cx, cy, cz, coarsecoords->getNumVectors());
308  this->writeFileVTKCells(cedgeStream, uniqueCoarseEdges, coarse_edges_vertices, coarse_edges_geomSize);
309  this->writeFileVTKClosing(cedgeStream);
310  cedgeStream.close();
311  }
312 #endif
313 
314  if (pL.get<int>("visualization: start level") <= fineLevel.GetLevelID()) {
315  // write out coarsening information
316  std::string filenameToWrite = this->getFileName(comm->getSize(), comm->getRank(), fineLevel.GetLevelID(), pL);
317  std::ofstream fout(filenameToWrite.c_str());
318 
319  std::vector<int> uniqueFine = this->makeUnique(vertices);
320  this->writeFileVTKOpening(fout, uniqueFine, geomSize);
321  this->writeFileVTKNodes(fout, uniqueFine, nodeMap);
322  this->writeFileVTKData(fout, uniqueFine, myAggOffset, vertex2AggId, comm->getRank());
323  this->writeFileVTKCoordinates(fout, uniqueFine, xCoords, yCoords, zCoords, coords->getNumVectors());
324  this->writeFileVTKCells(fout, uniqueFine, vertices, geomSize);
325  this->writeFileVTKClosing(fout);
326  fout.close();
327 
328  // create pvtu file
329  if (comm->getRank() == 0) {
330  std::string pvtuFilename = this->getPVTUFileName(comm->getSize(), comm->getRank(), fineLevel.GetLevelID(), pL);
331  std::string baseFname = this->getBaseFileName(comm->getSize(), fineLevel.GetLevelID(), pL);
332  std::ofstream pvtu(pvtuFilename.c_str());
333  this->writePVTU(pvtu, baseFname, comm->getSize(), pL.get<bool>("visualization: fine graph edges"));
334  pvtu.close();
335  }
336  }
337 
338  if (comm->getRank() == 0 && pL.get<bool>("visualization: build colormap")) {
339  this->buildColormap();
340  }
341 }
342 } // namespace MueLu
343 
344 #endif /* MUELU_AGGREGATIONEXPORTFACTORY_DEF_HPP_ */
Important warning messages (one line)
MueLu::DefaultLocalOrdinal LocalOrdinal
iterator begin() const
T & get(const std::string &name, T def_value)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
LocalOrdinal LO
ParameterList & set(std::string const &name, T &&value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
MueLu::DefaultNode Node
RCP< ParameterList > GetValidParameterList() const
bool isParameter(const std::string &name) const
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:63
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
const RCP< const Map > GetImportMap() const
Return overlapping import map (nodes).
iterator end() const
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)