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 // ***********************************************************************
4 //
5 // MueLu: A package for multigrid based preconditioning
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact
39 // Jonathan Hu (jhu@sandia.gov)
40 // Andrey Prokopenko (aprokop@sandia.gov)
41 // Ray Tuminaro (rstumin@sandia.gov)
42 //
43 // ***********************************************************************
44 //
45 // @HEADER
46 
47 #ifndef MUELU_COARSENINGVISUALIZATIONFACTORY_DEF_HPP_
48 #define MUELU_COARSENINGVISUALIZATIONFACTORY_DEF_HPP_
49 
50 #include <Xpetra_Matrix.hpp>
51 #include <Xpetra_ImportFactory.hpp>
52 #include <Xpetra_MultiVectorFactory.hpp>
53 #include "MueLu_CoarseningVisualizationFactory_decl.hpp"
54 #include "MueLu_Level.hpp"
55 
56 
57 namespace MueLu {
58 
59  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
60  RCP<const ParameterList> CoarseningVisualizationFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetValidParameterList() const {
61 
62  RCP<ParameterList> validParamList = VisualizationHelpers::GetValidParameterList();
63 
64  validParamList->set< int > ("visualization: start level", 0, "visualize only levels with level ids greater or equal than start level");// Remove me?
65 
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.");
70 
71  return validParamList;
72  }
73 
74  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
75  void CoarseningVisualizationFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::DeclareInput(Level &fineLevel, Level &coarseLevel) const {
76  this->Input(fineLevel, "Coordinates");
77 
78  const ParameterList & pL = this->GetParameterList();
79  TEUCHOS_TEST_FOR_EXCEPTION(pL.isParameter("P") && GetFactory("P") != Teuchos::null &&
80  pL.isParameter("Ptent") && GetFactory("Ptent") != Teuchos::null, Exceptions::RuntimeError,
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.");
84 
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");
89 
90  if(pL.get<bool>("visualization: fine graph edges"))
91  Input(fineLevel, "Graph");
92 #if 0
93  if(pL.get<bool>("visualization: coarse graph edges")) {
94  Input(coarseLevel, "Coordinates");
95  Input(coarseLevel, "Graph");
96  }
97 #endif
98  }
99 
100  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
101  void CoarseningVisualizationFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(Level &fineLevel, Level &coarseLevel) const {
102 
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");
110 
111  RCP<const Teuchos::Comm<int> > comm = P->getRowMap()->getComm();
112 
113  LocalOrdinal dofsPerNode = 1;
114  LocalOrdinal stridedRowOffset = 0;
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"));
118  LocalOrdinal blockid = strRowMap->getStridedBlockId();
119  if (blockid > -1) {
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]);
124  } else {
125  dofsPerNode = strRowMap->getFixedBlockSize();
126  }
127  GetOStream(Runtime1) << "CoarseningVisualizationFactory::Build():" << " #dofs per node = " << dofsPerNode << std::endl;
128  }
129 
130  LocalOrdinal columnsPerNode = dofsPerNode;
131  LocalOrdinal stridedColumnOffset = 0;
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();
136 
137  if (blockid > -1) {
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]);
142  } else {
143  columnsPerNode = strDomainMap->getFixedBlockSize();
144  }
145  GetOStream(Runtime1) << "CoarseningVisualizationFactory::Build():" << " #columns per node = " << columnsPerNode << std::endl;
146  }
147 
148  // TODO add support for overlapping aggregates
149  TEUCHOS_TEST_FOR_EXCEPTION(strDomainMap->getNodeNumElements() != P->getColMap()->getNodeNumElements(), Exceptions::RuntimeError,
150  "CoarseningVisualization only supports non-overlapping transfers");
151 
152  // number of local "aggregates"
153  LocalOrdinal numLocalAggs = strDomainMap->getNodeNumElements() / columnsPerNode;
154  std::vector< std::set<LocalOrdinal> > localAggs(numLocalAggs);
155 
156  // do loop over all local rows of prolongator and extract column information
157  for (LO row = 0; row < Teuchos::as<LO>(P->getRowMap()->getNodeNumElements()); ++row) {
158  ArrayView<const LO> indices;
159  ArrayView<const SC> vals;
160  P->getLocalRowView(row, indices, vals);
161 
162  for(typename ArrayView<const LO>::iterator c = indices.begin(); c != indices.end(); ++c) {
163  localAggs[(*c)/columnsPerNode].insert(row/dofsPerNode);
164  }
165  }
166 
167  // determine number of "aggs" per proc and calculate local "agg" offset...
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]);
172 
173  LocalOrdinal myAggOffset = 0;
174  for(int i = 0; i < comm->getRank(); ++i) {
175  myAggOffset += numLocalAggsPerProc[i];
176  }
177 
178  /*for (LocalOrdinal i = 0; i < numLocalAggs; ++i) {
179 
180  std::cout << "PROC: " << comm->getRank() << " Local aggregate: " << i + myAggOffset << " with nodes: ";
181  for( typename std::set<LocalOrdinal>::iterator it = localAggs[i].begin(); it != localAggs[i].end(); ++it) {
182  std::cout << *it << ", ";
183  }
184  std::cout << std::endl;
185  }*/
186 
187  // get fine level coordinate information
188  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");
189 
190  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::as<LO>(P->getRowMap()->getNodeNumElements()) / dofsPerNode != Teuchos::as<LocalOrdinal>(coords->getLocalLength()), Exceptions::RuntimeError,
191  "Number of fine level nodes in coordinates is inconsistent with dof based information");
192 
193  // communicate fine level coordinates if necessary
194  if (pL.get<bool>("visualization: fine graph edges")) {
195  fineGraph = Get<RCP<GraphBase> >(fineLevel, "Graph");
196 
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;
201  }
202 
203  Teuchos::RCP<const Map> nodeMap = coords->getMap();
204 
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));
210  }
211 
212  // determine number of nodes on fine level
213  LocalOrdinal numFineNodes = Teuchos::as<LocalOrdinal>(coords->getLocalLength());
214 
215  // create vertex2AggId array
216  ArrayRCP<LocalOrdinal> vertex2AggId(numFineNodes, -1);
217  for (LocalOrdinal i = 0; i < numLocalAggs; ++i) {
218  // TODO: check if entry = -1
219  for( typename std::set<LocalOrdinal>::iterator it = localAggs[i].begin(); it != localAggs[i].end(); ++it) {
220  vertex2AggId[*it] = i;
221  }
222  }
223 
224  // we have no information which node is root and which is not
225  // we could either look at the entries in P again or build some new root nodes
226  // assuming that root nodes are usually at the center of the aggregate
227  std::vector<bool> isRoot(numFineNodes, false);
228  for (LocalOrdinal i = 0; i < numLocalAggs; ++i) {
229 
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;
233 
234  // loop over all nodes in aggregate i and determine center coordinates of aggregate
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];
239  }
240  xCenter /= Teuchos::as<LocalOrdinal>(localAggs[i].size());
241  yCenter /= Teuchos::as<LocalOrdinal>(localAggs[i].size());
242  zCenter /= Teuchos::as<LocalOrdinal>(localAggs[i].size());
243 
244  // loop over all nodes in aggregate i and find node which is closest to aggregate center
245  LocalOrdinal rootCandidate = -1;
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;
259  rootCandidate = *it;
260  }
261  }
262 
263  isRoot[rootCandidate] = true;
264  }
265 
266  std::vector<LocalOrdinal> vertices;
267  std::vector<LocalOrdinal> geomSize;
268  int vizLevel = pL.get<int>("visualization: start level");
269  if(vizLevel <= fineLevel.GetLevelID()) {
270 
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") {
277  // TODO do a smarter distinction and check the number of z levels...
278  // loop over all coordinates and collect coordinate components in sets...
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);
283  }
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);
290 #endif
291  }
292  else
293  {
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);
297  }
298  }
299 
300  // write out fine edge information
301  if(pL.get<bool>("visualization: fine graph edges")) {
302  TEUCHOS_TEST_FOR_EXCEPTION(fineGraph == Teuchos::null, Exceptions::RuntimeError,
303  "Could not get information about fine graph.");
304 
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);
308 
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());
312 
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);
320  edgeStream.close();
321  }
322 
323  // communicate fine level coordinates if necessary
324 #if 0 // we don't have access to the coarse graph
325  if (pL.get<bool>("visualization: coarse graph edges")) {
326  RCP<GraphBase> coarseGraph = Get<RCP<GraphBase> >(coarseLevel, "Graph");
327 
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");
329 
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;
334 
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));
340  }
341 
342  Teuchos::RCP<const Map> coarsenodeMap = coarsecoords->getMap();
343 
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);
347 
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());
351 
352  std::vector<int> uniqueCoarseEdges = this->makeUnique(coarse_edges_vertices);
353  this->writeFileVTKOpening(cedgeStream, uniqueCoarseEdges, coarse_edges_geomSize);
354  this->writeFileVTKNodes(cedgeStream, uniqueCoarseEdges, coarsenodeMap);
355  //this->writeFileVTKData(edgeStream, uniqueCoarseEdges, myAggOffset, vertex2AggId, comm->getRank());
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);
359  cedgeStream.close();
360  }
361 #endif
362 
363  if(pL.get<int>("visualization: start level") <= fineLevel.GetLevelID()) {
364  // write out coarsening information
365  std::string filenameToWrite = this->getFileName(comm->getSize(), comm->getRank(), fineLevel.GetLevelID(), pL);
366  std::ofstream fout (filenameToWrite.c_str());
367 
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);
375  fout.close();
376 
377  // create pvtu file
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"));
383  pvtu.close();
384  }
385  }
386 
387  if(comm->getRank() == 0 && pL.get<bool>("visualization: build colormap")) {
388  this->buildColormap();
389  }
390  }
391 } // namespace MueLu
392 
393 #endif /* MUELU_AGGREGATIONEXPORTFACTORY_DEF_HPP_ */
MueLu::Warnings0
Important warning messages (one line)
Definition: MueLu_VerbosityLevel.hpp:57
Xpetra::ImportFactory::Build
static RCP< Import< LocalOrdinal, GlobalOrdinal, Node > > Build(const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &source, const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &target)
Teuchos::ArrayView::begin
iterator begin() const
Teuchos::ParameterList::get
T & get(const std::string &name, T def_value)
Teuchos::ParameterList::set
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
LO
LocalOrdinal LO
MueLu_Level.hpp
MueLu::CoarseningVisualizationFactory::GetValidParameterList
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Definition: MueLu_CoarseningVisualizationFactory_def.hpp:60
Xpetra_Matrix.hpp
MueLu::VisualizationHelpers::GetValidParameterList
RCP< ParameterList > GetValidParameterList() const
Definition: MueLu_VisualizationHelpers_def.hpp:74
Teuchos::ParameterList::isParameter
bool isParameter(const std::string &name) const
Xpetra_MultiVectorFactory.hpp
MueLu::Level
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
Teuchos::reduceAll
TEUCHOS_DEPRECATED void reduceAll(const Comm< Ordinal > &comm, const EReductionType reductType, const Packet &send, Packet *globalReduct)
MueLu_CoarseningVisualizationFactory_decl.hpp
MueLu::CoarseningVisualizationFactory::DeclareInput
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
Definition: MueLu_CoarseningVisualizationFactory_def.hpp:75
Teuchos::ArrayView::end
iterator end() const
MueLu::GraphBase::GetImportMap
virtual const RCP< const Map > GetImportMap() const =0
MueLu::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition: MueLu_Exceptions.hpp:70
MueLu::Runtime1
Description of what is happening (more verbose)
Definition: MueLu_VerbosityLevel.hpp:63
MueLu::CoarseningVisualizationFactory::Build
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
Definition: MueLu_CoarseningVisualizationFactory_def.hpp:101
Xpetra_ImportFactory.hpp