MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_AggregationExportFactory_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 /*
11  * MueLu_AggregationExportFactory_def.hpp
12  *
13  * Created on: Feb 10, 2012
14  * Author: wiesner
15  */
16 
17 #ifndef MUELU_AGGREGATIONEXPORTFACTORY_DEF_HPP_
18 #define MUELU_AGGREGATIONEXPORTFACTORY_DEF_HPP_
19 
20 #include <Xpetra_Matrix.hpp>
21 #include <Xpetra_CrsMatrixWrap.hpp>
22 #include <Xpetra_ImportFactory.hpp>
23 #include <Xpetra_MultiVectorFactory.hpp>
24 #include "MueLu_AggregationExportFactory_decl.hpp"
25 #include "MueLu_Level.hpp"
26 #include "MueLu_Aggregates.hpp"
27 
28 #include "MueLu_AmalgamationFactory.hpp"
29 #include "MueLu_AmalgamationInfo.hpp"
30 #include "MueLu_Monitor.hpp"
31 #include <vector>
32 #include <list>
33 #include <algorithm>
34 #include <string>
35 #include <stdexcept>
36 #include <cstdio>
37 #include <cmath>
38 
39 namespace MueLu {
40 
41 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
42 AggregationExportFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::AggregationExportFactory()
43  : doFineGraphEdges_(false)
44  , doCoarseGraphEdges_(false)
45  , numNodes_(0)
46  , numAggs_(0)
47  , dims_(0)
48  , myRank_(-1)
49  , aggsOffset_(0) {}
50 
51 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
52 AggregationExportFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::~AggregationExportFactory() = default;
53 
54 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
55 RCP<const ParameterList> AggregationExportFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetValidParameterList() const {
56  RCP<ParameterList> validParamList = rcp(new ParameterList());
57 
58  std::string output_msg =
59  "Output filename template (%TIMESTEP is replaced by \'Output file: time step\' variable,"
60  "%ITER is replaced by \'Output file: iter\' variable, %LEVELID is replaced level id, %PROCID is replaced by processor id)";
61  std::string output_def = "aggs_level%LEVELID_proc%PROCID.out";
62 
63  validParamList->set<RCP<const FactoryBase> >("A", Teuchos::null, "Factory for A.");
64  validParamList->set<RCP<const FactoryBase> >("Coordinates", Teuchos::null, "Factory for Coordinates.");
65  validParamList->set<RCP<const FactoryBase> >("Graph", Teuchos::null, "Factory for Graph.");
66  validParamList->set<RCP<const FactoryBase> >("Aggregates", Teuchos::null, "Factory for Aggregates.");
67  validParamList->set<RCP<const FactoryBase> >("UnAmalgamationInfo", Teuchos::null, "Factory for UnAmalgamationInfo.");
68  validParamList->set<RCP<const FactoryBase> >("DofsPerNode", Teuchos::null, "Factory for DofsPerNode.");
69  // CMS/BMK: Old style factory-only options. Deprecate me.
70  validParamList->set<std::string>("Output filename", output_def, output_msg);
71  validParamList->set<int>("Output file: time step", 0, "time step variable for output file name");
72  validParamList->set<int>("Output file: iter", 0, "nonlinear iteration variable for output file name");
73 
74  // New-style master list options (here are same defaults as in masterList.xml)
75  validParamList->set<std::string>("aggregation: output filename", "", "filename for VTK-style visualization output");
76  validParamList->set<int>("aggregation: output file: time step", 0, "time step variable for output file name"); // Remove me?
77  validParamList->set<int>("aggregation: output file: iter", 0, "nonlinear iteration variable for output file name"); // Remove me?
78  validParamList->set<std::string>("aggregation: output file: agg style", "Point Cloud", "style of aggregate visualization for VTK output");
79  validParamList->set<bool>("aggregation: output file: fine graph edges", false, "Whether to draw all fine node connections along with the aggregates.");
80  validParamList->set<bool>("aggregation: output file: coarse graph edges", false, "Whether to draw all coarse node connections along with the aggregates.");
81  validParamList->set<bool>("aggregation: output file: build colormap", false, "Whether to output a random colormap for ParaView in a separate XML file.");
82  return validParamList;
83 }
84 
85 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
86 void AggregationExportFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::DeclareInput(Level& fineLevel, Level& coarseLevel) const {
87  Input(fineLevel, "Aggregates"); //< factory which created aggregates
88  Input(fineLevel, "DofsPerNode"); //< CoalesceAndDropFactory (needed for DofsPerNode variable)
89  Input(fineLevel, "UnAmalgamationInfo"); //< AmalgamationFactory (needed for UnAmalgamationInfo variable)
90 
91  const ParameterList& pL = GetParameterList();
92  // Only pull in coordinates if the user explicitly requests direct VTK output, so as not to break uses of old code
93  if (pL.isParameter("aggregation: output filename") && pL.get<std::string>("aggregation: output filename").length()) {
94  Input(fineLevel, "Coordinates");
95  Input(fineLevel, "A");
96  Input(fineLevel, "Graph");
97  if (pL.get<bool>("aggregation: output file: coarse graph edges")) {
98  Input(coarseLevel, "Coordinates");
99  Input(coarseLevel, "A");
100  Input(coarseLevel, "Graph");
101  }
102  }
103 }
104 
105 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
106 void AggregationExportFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(Level& fineLevel, Level& coarseLevel) const {
107  using namespace std;
108  // Decide which build function to follow, based on input params
109  const ParameterList& pL = GetParameterList();
110  FactoryMonitor m(*this, "AggregationExportFactory", coarseLevel);
111  Teuchos::RCP<Aggregates> aggregates = Get<Teuchos::RCP<Aggregates> >(fineLevel, "Aggregates");
112  Teuchos::RCP<const Teuchos::Comm<int> > comm = aggregates->GetMap()->getComm();
113  int numProcs = comm->getSize();
114  int myRank = comm->getRank();
115  string masterFilename = pL.get<std::string>("aggregation: output filename"); // filename parameter from master list
116  string pvtuFilename = ""; // only root processor will set this
117  string localFilename = pL.get<std::string>("Output filename");
118  string filenameToWrite;
119  bool useVTK = false;
120  doCoarseGraphEdges_ = pL.get<bool>("aggregation: output file: coarse graph edges");
121  doFineGraphEdges_ = pL.get<bool>("aggregation: output file: fine graph edges");
122  if (masterFilename.length()) {
123  useVTK = true;
124  filenameToWrite = masterFilename;
125  if (filenameToWrite.rfind(".vtu") == string::npos) // Must have the file extension in the name
126  filenameToWrite.append(".vtu");
127  if (numProcs > 1 && filenameToWrite.rfind("%PROCID") == string::npos) // filename can't be identical between processsors in parallel problem
128  filenameToWrite.insert(filenameToWrite.rfind(".vtu"), "-proc%PROCID");
129  } else
130  filenameToWrite = localFilename;
131  LocalOrdinal DofsPerNode = Get<LocalOrdinal>(fineLevel, "DofsPerNode");
132  Teuchos::RCP<AmalgamationInfo> amalgInfo = Get<RCP<AmalgamationInfo> >(fineLevel, "UnAmalgamationInfo");
133  Teuchos::RCP<Matrix> Amat = Get<RCP<Matrix> >(fineLevel, "A");
134  Teuchos::RCP<Matrix> Ac;
135  if (doCoarseGraphEdges_)
136  Ac = Get<RCP<Matrix> >(coarseLevel, "A");
137  Teuchos::RCP<CoordinateMultiVector> coords = Teuchos::null;
138  Teuchos::RCP<CoordinateMultiVector> coordsCoarse = Teuchos::null;
139  Teuchos::RCP<LWGraph> fineGraph = Teuchos::null;
140  Teuchos::RCP<LWGraph> coarseGraph = Teuchos::null;
141  if (doFineGraphEdges_)
142  fineGraph = Get<RCP<LWGraph> >(fineLevel, "Graph");
143  if (doCoarseGraphEdges_)
144  coarseGraph = Get<RCP<LWGraph> >(coarseLevel, "Graph");
145  if (useVTK) // otherwise leave null, will not be accessed by non-vtk code
146  {
147  coords = Get<RCP<CoordinateMultiVector> >(fineLevel, "Coordinates");
148  coords_ = coords;
149  if (doCoarseGraphEdges_)
150  coordsCoarse = Get<RCP<CoordinateMultiVector> >(coarseLevel, "Coordinates");
151  dims_ = coords->getNumVectors(); // 2D or 3D?
152  if (numProcs > 1) {
153  if (aggregates->AggregatesCrossProcessors()) { // Do we want to use the map from aggregates here instead of the map from A? Using the map from A seems to be problematic with multiple dofs per node
154  RCP<Import> coordImporter = Xpetra::ImportFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(coords->getMap(), Amat->getColMap());
155  RCP<CoordinateMultiVector> ghostedCoords = Xpetra::MultiVectorFactory<coordinate_type, LocalOrdinal, GlobalOrdinal, Node>::Build(Amat->getColMap(), dims_);
156  ghostedCoords->doImport(*coords, *coordImporter, Xpetra::INSERT);
157  coords = ghostedCoords;
158  coords_ = ghostedCoords;
159  }
160  if (doCoarseGraphEdges_) {
161  RCP<Import> coordImporter = Xpetra::ImportFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(coordsCoarse->getMap(), Ac->getColMap());
162  RCP<CoordinateMultiVector> ghostedCoords = Xpetra::MultiVectorFactory<coordinate_type, LocalOrdinal, GlobalOrdinal, Node>::Build(Ac->getColMap(), dims_);
163  ghostedCoords->doImport(*coordsCoarse, *coordImporter, Xpetra::INSERT);
164  coordsCoarse = ghostedCoords;
165  coordsCoarse_ = ghostedCoords;
166  }
167  }
168  }
169  GetOStream(Runtime0) << "AggregationExportFactory: DofsPerNode: " << DofsPerNode << std::endl;
170  Teuchos::RCP<LocalOrdinalMultiVector> vertex2AggId_vector = aggregates->GetVertex2AggId();
171  Teuchos::RCP<LocalOrdinalVector> procWinner_vector = aggregates->GetProcWinner();
172  Teuchos::ArrayRCP<LocalOrdinal> vertex2AggId = aggregates->GetVertex2AggId()->getDataNonConst(0);
173  Teuchos::ArrayRCP<LocalOrdinal> procWinner = aggregates->GetProcWinner()->getDataNonConst(0);
174 
175  vertex2AggId_ = vertex2AggId;
176 
177  // prepare for calculating global aggregate ids
178  std::vector<GlobalOrdinal> numAggsGlobal(numProcs, 0);
179  std::vector<GlobalOrdinal> numAggsLocal(numProcs, 0);
180  std::vector<GlobalOrdinal> minGlobalAggId(numProcs, 0);
181 
182  numAggsLocal[myRank] = aggregates->GetNumAggregates();
183  Teuchos::reduceAll(*comm, Teuchos::REDUCE_SUM, numProcs, &numAggsLocal[0], &numAggsGlobal[0]);
184  for (int i = 1; i < Teuchos::as<int>(numAggsGlobal.size()); ++i) {
185  numAggsGlobal[i] += numAggsGlobal[i - 1];
186  minGlobalAggId[i] = numAggsGlobal[i - 1];
187  }
188  if (numProcs == 0)
189  aggsOffset_ = 0;
190  else
191  aggsOffset_ = minGlobalAggId[myRank];
192  ArrayRCP<LO> aggStart;
193  ArrayRCP<GlobalOrdinal> aggToRowMap;
194  amalgInfo->UnamalgamateAggregates(*aggregates, aggStart, aggToRowMap);
195  int timeStep = pL.get<int>("Output file: time step");
196  int iter = pL.get<int>("Output file: iter");
197  filenameToWrite = this->replaceAll(filenameToWrite, "%LEVELID", toString(fineLevel.GetLevelID()));
198  filenameToWrite = this->replaceAll(filenameToWrite, "%TIMESTEP", toString(timeStep));
199  filenameToWrite = this->replaceAll(filenameToWrite, "%ITER", toString(iter));
200  // Proc id MUST be included in vtu filenames to distinguish them (if multiple procs)
201  // In all other cases (else), including processor # in filename is optional
202  string masterStem = "";
203  if (useVTK) {
204  masterStem = filenameToWrite.substr(0, filenameToWrite.rfind(".vtu"));
205  masterStem = this->replaceAll(masterStem, "%PROCID", "");
206  }
207  pvtuFilename = masterStem + "-master.pvtu";
208  string baseFname = filenameToWrite; // get a version of the filename string with the %PROCID token, but without substituting myRank (needed for pvtu output)
209  filenameToWrite = this->replaceAll(filenameToWrite, "%PROCID", toString(myRank));
210  GetOStream(Runtime0) << "AggregationExportFactory: outputfile \"" << filenameToWrite << "\"" << std::endl;
211  ofstream fout(filenameToWrite.c_str());
212  GO numAggs = aggregates->GetNumAggregates();
213  if (!useVTK) {
214  GO indexBase = aggregates->GetMap()->getIndexBase(); // extract indexBase from overlapping map within aggregates structure. The indexBase is constant throughout the whole simulation (either 0 = C++ or 1 = Fortran)
215  GO offset = amalgInfo->GlobalOffset(); // extract offset for global dof ids
216  vector<GlobalOrdinal> nodeIds;
217  for (int i = 0; i < numAggs; ++i) {
218  fout << "Agg " << minGlobalAggId[myRank] + i << " Proc " << myRank << ":";
219 
220  // TODO: Use k+=DofsPerNode instead of ++k and get rid of std::unique call afterwards
221  for (int k = aggStart[i]; k < aggStart[i + 1]; ++k) {
222  nodeIds.push_back((aggToRowMap[k] - offset - indexBase) / DofsPerNode + indexBase);
223  }
224 
225  // remove duplicate entries from nodeids
226  std::sort(nodeIds.begin(), nodeIds.end());
227  typename std::vector<GlobalOrdinal>::iterator endLocation = std::unique(nodeIds.begin(), nodeIds.end());
228  nodeIds.erase(endLocation, nodeIds.end());
229 
230  // print out nodeids
231  for (typename std::vector<GlobalOrdinal>::iterator printIt = nodeIds.begin(); printIt != nodeIds.end(); printIt++)
232  fout << " " << *printIt;
233  nodeIds.clear();
234  fout << std::endl;
235  }
236  } else {
237  using namespace std;
238  // Make sure we have coordinates
239  TEUCHOS_TEST_FOR_EXCEPTION(coords.is_null(), Exceptions::RuntimeError, "AggExportFactory could not get coordinates, but they are required for VTK output.");
240  numAggs_ = numAggs;
241  numNodes_ = coords->getLocalLength();
242  // Get the sizes of the aggregates to speed up grabbing node IDs
243  aggSizes_ = aggregates->ComputeAggregateSizesArrayRCP();
244  myRank_ = myRank;
245  string aggStyle = "Point Cloud";
246  try {
247  aggStyle = pL.get<string>("aggregation: output file: agg style"); // Let "Point Cloud" be the default style
248  } catch (std::exception& e) {
249  }
250  vector<int> vertices;
251  vector<int> geomSizes;
252  string indent = "";
253  nodeMap_ = Amat->getMap();
254  for (LocalOrdinal i = 0; i < numNodes_; i++) {
255  isRoot_.push_back(aggregates->IsRoot(i));
256  }
257  // If problem is serial and not outputting fine nor coarse graph edges, don't make pvtu file
258  // Otherwise create it
259  if (myRank == 0 && (numProcs != 1 || doCoarseGraphEdges_ || doFineGraphEdges_)) {
260  ofstream pvtu(pvtuFilename.c_str());
261  writePVTU_(pvtu, baseFname, numProcs);
262  pvtu.close();
263  }
264  if (aggStyle == "Point Cloud")
265  this->doPointCloud(vertices, geomSizes, numAggs_, numNodes_);
266  else if (aggStyle == "Jacks")
267  this->doJacks(vertices, geomSizes, numAggs_, numNodes_, isRoot_, vertex2AggId_);
268  else if (aggStyle == "Jacks++") // Not actually implemented
269  doJacksPlus_(vertices, geomSizes);
270  else if (aggStyle == "Convex Hulls")
271  doConvexHulls(vertices, geomSizes);
272  else {
273  GetOStream(Warnings0) << " Unrecognized agg style.\n Possible values are Point Cloud, Jacks, Jacks++, and Convex Hulls.\n Defaulting to Point Cloud." << std::endl;
274  aggStyle = "Point Cloud";
275  this->doPointCloud(vertices, geomSizes, numAggs_, numNodes_);
276  }
277  writeFile_(fout, aggStyle, vertices, geomSizes);
278  if (doCoarseGraphEdges_) {
279  string fname = filenameToWrite;
280  string cEdgeFile = fname.insert(fname.rfind(".vtu"), "-coarsegraph");
281  std::ofstream edgeStream(cEdgeFile.c_str());
282  doGraphEdges_(edgeStream, Ac, coarseGraph, false, DofsPerNode);
283  }
284  if (doFineGraphEdges_) {
285  string fname = filenameToWrite;
286  string fEdgeFile = fname.insert(fname.rfind(".vtu"), "-finegraph");
287  std::ofstream edgeStream(fEdgeFile.c_str());
288  doGraphEdges_(edgeStream, Amat, fineGraph, true, DofsPerNode);
289  }
290  if (myRank == 0 && pL.get<bool>("aggregation: output file: build colormap")) {
291  buildColormap_();
292  }
293  }
294  fout.close();
295 }
296 
297 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
298 void AggregationExportFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::doJacksPlus_(std::vector<int>& /* vertices */, std::vector<int>& /* geomSizes */) const {
299  // TODO
300 }
301 
302 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
303 void AggregationExportFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::doConvexHulls(std::vector<int>& vertices, std::vector<int>& geomSizes) const {
304  Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> xCoords = coords_->getData(0);
305  Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> yCoords = coords_->getData(1);
306  Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> zCoords = Teuchos::null;
307 
308  if (dims_ == 2) {
309  this->doConvexHulls2D(vertices, geomSizes, numAggs_, numNodes_, isRoot_, vertex2AggId_, xCoords, yCoords);
310  } else {
311  zCoords = coords_->getData(2);
312  this->doConvexHulls3D(vertices, geomSizes, numAggs_, numNodes_, isRoot_, vertex2AggId_, xCoords, yCoords, zCoords);
313  }
314 }
315 
316 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
317 void AggregationExportFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::doGraphEdges_(std::ofstream& fout, Teuchos::RCP<Matrix>& A, Teuchos::RCP<LWGraph>& G, bool fine, int dofs) const {
318  using namespace std;
319  ArrayView<const Scalar> values;
320  // Allow two different colors of connections (by setting "aggregates" scalar to CONTRAST_1 or CONTRAST_2)
321  vector<pair<int, int> > vert1; // vertices (node indices)
322  vector<pair<int, int> > vert2; // size of every cell is assumed to be 2 vertices, since all edges are drawn as lines
323 
324  Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> xCoords = coords_->getData(0);
325  Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> yCoords = coords_->getData(1);
326  Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> zCoords = Teuchos::null;
327  if (dims_ == 3)
328  zCoords = coords_->getData(2);
329 
330  Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> cx = coordsCoarse_->getData(0);
331  Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> cy = coordsCoarse_->getData(1);
332  Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> cz = Teuchos::null;
333  if (dims_ == 3)
334  cz = coordsCoarse_->getData(2);
335 
336  if (A->isGloballyIndexed()) {
337  ArrayView<const GlobalOrdinal> indices;
338  for (GlobalOrdinal globRow = 0; globRow < GlobalOrdinal(A->getGlobalNumRows()); globRow++) {
339  if (dofs == 1)
340  A->getGlobalRowView(globRow, indices, values);
341  auto neighbors = G->getNeighborVertices((LocalOrdinal)globRow);
342  int gEdge = 0;
343  int aEdge = 0;
344  while (gEdge != int(neighbors.length)) {
345  if (dofs == 1) {
346  if (neighbors(gEdge) == indices[aEdge]) {
347  // graph and matrix both have this edge, wasn't filtered, show as color 1
348  vert1.push_back(pair<int, int>(int(globRow), neighbors(gEdge)));
349  gEdge++;
350  aEdge++;
351  } else {
352  // graph contains an edge at gEdge which was filtered from A, show as color 2
353  vert2.push_back(pair<int, int>(int(globRow), neighbors(gEdge)));
354  gEdge++;
355  }
356  } else // for multiple DOF problems, don't try to detect filtered edges and ignore A
357  {
358  vert1.push_back(pair<int, int>(int(globRow), neighbors(gEdge)));
359  gEdge++;
360  }
361  }
362  }
363  } else {
364  ArrayView<const LocalOrdinal> indices;
365  for (LocalOrdinal locRow = 0; locRow < LocalOrdinal(A->getLocalNumRows()); locRow++) {
366  if (dofs == 1)
367  A->getLocalRowView(locRow, indices, values);
368  auto neighbors = G->getNeighborVertices(locRow);
369  // Add those local indices (columns) to the list of connections (which are pairs of the form (localM, localN))
370  int gEdge = 0;
371  int aEdge = 0;
372  while (gEdge != int(neighbors.length)) {
373  if (dofs == 1) {
374  if (neighbors(gEdge) == indices[aEdge]) {
375  vert1.push_back(pair<int, int>(locRow, neighbors(gEdge)));
376  gEdge++;
377  aEdge++;
378  } else {
379  vert2.push_back(pair<int, int>(locRow, neighbors(gEdge)));
380  gEdge++;
381  }
382  } else {
383  vert1.push_back(pair<int, int>(locRow, neighbors(gEdge)));
384  gEdge++;
385  }
386  }
387  }
388  }
389  for (size_t i = 0; i < vert1.size(); i++) {
390  if (vert1[i].first > vert1[i].second) {
391  int temp = vert1[i].first;
392  vert1[i].first = vert1[i].second;
393  vert1[i].second = temp;
394  }
395  }
396  for (size_t i = 0; i < vert2.size(); i++) {
397  if (vert2[i].first > vert2[i].second) {
398  int temp = vert2[i].first;
399  vert2[i].first = vert2[i].second;
400  vert2[i].second = temp;
401  }
402  }
403  sort(vert1.begin(), vert1.end());
404  vector<pair<int, int> >::iterator newEnd = unique(vert1.begin(), vert1.end()); // remove duplicate edges
405  vert1.erase(newEnd, vert1.end());
406  sort(vert2.begin(), vert2.end());
407  newEnd = unique(vert2.begin(), vert2.end());
408  vert2.erase(newEnd, vert2.end());
409  vector<int> points1;
410  points1.reserve(2 * vert1.size());
411  for (size_t i = 0; i < vert1.size(); i++) {
412  points1.push_back(vert1[i].first);
413  points1.push_back(vert1[i].second);
414  }
415  vector<int> points2;
416  points2.reserve(2 * vert2.size());
417  for (size_t i = 0; i < vert2.size(); i++) {
418  points2.push_back(vert2[i].first);
419  points2.push_back(vert2[i].second);
420  }
421  vector<int> unique1 = this->makeUnique(points1);
422  vector<int> unique2 = this->makeUnique(points2);
423  fout << "<VTKFile type=\"UnstructuredGrid\" byte_order=\"LittleEndian\">" << endl;
424  fout << " <UnstructuredGrid>" << endl;
425  fout << " <Piece NumberOfPoints=\"" << unique1.size() + unique2.size() << "\" NumberOfCells=\"" << vert1.size() + vert2.size() << "\">" << endl;
426  fout << " <PointData Scalars=\"Node Aggregate Processor\">" << endl;
427  fout << " <DataArray type=\"Int32\" Name=\"Node\" format=\"ascii\">" << endl; // node and aggregate will be set to CONTRAST_1|2, but processor will have its actual value
428  string indent = " ";
429  fout << indent;
430  for (size_t i = 0; i < unique1.size(); i++) {
431  fout << CONTRAST_1_ << " ";
432  if (i % 25 == 24)
433  fout << endl
434  << indent;
435  }
436  for (size_t i = 0; i < unique2.size(); i++) {
437  fout << CONTRAST_2_ << " ";
438  if ((i + 2 * vert1.size()) % 25 == 24)
439  fout << endl
440  << indent;
441  }
442  fout << endl;
443  fout << " </DataArray>" << endl;
444  fout << " <DataArray type=\"Int32\" Name=\"Aggregate\" format=\"ascii\">" << endl;
445  fout << indent;
446  for (size_t i = 0; i < unique1.size(); i++) {
447  fout << CONTRAST_1_ << " ";
448  if (i % 25 == 24)
449  fout << endl
450  << indent;
451  }
452  for (size_t i = 0; i < unique2.size(); i++) {
453  fout << CONTRAST_2_ << " ";
454  if ((i + 2 * vert2.size()) % 25 == 24)
455  fout << endl
456  << indent;
457  }
458  fout << endl;
459  fout << " </DataArray>" << endl;
460  fout << " <DataArray type=\"Int32\" Name=\"Processor\" format=\"ascii\">" << endl;
461  fout << indent;
462  for (size_t i = 0; i < unique1.size() + unique2.size(); i++) {
463  fout << myRank_ << " ";
464  if (i % 25 == 24)
465  fout << endl
466  << indent;
467  }
468  fout << endl;
469  fout << " </DataArray>" << endl;
470  fout << " </PointData>" << endl;
471  fout << " <Points>" << endl;
472  fout << " <DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"ascii\">" << endl;
473  fout << indent;
474  for (size_t i = 0; i < unique1.size(); i++) {
475  if (fine) {
476  fout << xCoords[unique1[i]] << " " << yCoords[unique1[i]] << " ";
477  if (dims_ == 3)
478  fout << zCoords[unique1[i]] << " ";
479  else
480  fout << "0 ";
481  if (i % 2)
482  fout << endl
483  << indent;
484  } else {
485  fout << cx[unique1[i]] << " " << cy[unique1[i]] << " ";
486  if (dims_ == 3)
487  fout << cz[unique1[i]] << " ";
488  else
489  fout << "0 ";
490  if (i % 2)
491  fout << endl
492  << indent;
493  }
494  }
495  for (size_t i = 0; i < unique2.size(); i++) {
496  if (fine) {
497  fout << xCoords[unique2[i]] << " " << yCoords[unique2[i]] << " ";
498  if (dims_ == 3)
499  fout << zCoords[unique2[i]] << " ";
500  else
501  fout << "0 ";
502  if (i % 2)
503  fout << endl
504  << indent;
505  } else {
506  fout << cx[unique2[i]] << " " << cy[unique2[i]] << " ";
507  if (dims_ == 3)
508  fout << cz[unique2[i]] << " ";
509  else
510  fout << "0 ";
511  if ((i + unique1.size()) % 2)
512  fout << endl
513  << indent;
514  }
515  }
516  fout << endl;
517  fout << " </DataArray>" << endl;
518  fout << " </Points>" << endl;
519  fout << " <Cells>" << endl;
520  fout << " <DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">" << endl;
521  fout << indent;
522  for (size_t i = 0; i < points1.size(); i++) {
523  fout << points1[i] << " ";
524  if (i % 10 == 9)
525  fout << endl
526  << indent;
527  }
528  for (size_t i = 0; i < points2.size(); i++) {
529  fout << points2[i] + unique1.size() << " ";
530  if ((i + points1.size()) % 10 == 9)
531  fout << endl
532  << indent;
533  }
534  fout << endl;
535  fout << " </DataArray>" << endl;
536  fout << " <DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">" << endl;
537  fout << indent;
538  int offset = 0;
539  for (size_t i = 0; i < vert1.size() + vert2.size(); i++) {
540  offset += 2;
541  fout << offset << " ";
542  if (i % 10 == 9)
543  fout << endl
544  << indent;
545  }
546  fout << endl;
547  fout << " </DataArray>" << endl;
548  fout << " <DataArray type=\"Int32\" Name=\"types\" format=\"ascii\">" << endl;
549  fout << indent;
550  for (size_t i = 0; i < vert1.size() + vert2.size(); i++) {
551  fout << "3 ";
552  if (i % 25 == 24)
553  fout << endl
554  << indent;
555  }
556  fout << endl;
557  fout << " </DataArray>" << endl;
558  fout << " </Cells>" << endl;
559  fout << " </Piece>" << endl;
560  fout << " </UnstructuredGrid>" << endl;
561  fout << "</VTKFile>" << endl;
562 }
563 
564 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
565 void AggregationExportFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeFile_(std::ofstream& fout, std::string styleName, std::vector<int>& vertices, std::vector<int>& geomSizes) const {
566  using namespace std;
567 
568  Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> xCoords = coords_->getData(0);
569  Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> yCoords = coords_->getData(1);
570  Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> zCoords = Teuchos::null;
571  if (dims_ == 3)
572  zCoords = coords_->getData(2);
573 
574  vector<int> uniqueFine = this->makeUnique(vertices);
575  string indent = " ";
576  fout << "<!--" << styleName << " Aggregates Visualization-->" << endl;
577  fout << "<VTKFile type=\"UnstructuredGrid\" byte_order=\"LittleEndian\">" << endl;
578  fout << " <UnstructuredGrid>" << endl;
579  fout << " <Piece NumberOfPoints=\"" << uniqueFine.size() << "\" NumberOfCells=\"" << geomSizes.size() << "\">" << endl;
580  fout << " <PointData Scalars=\"Node Aggregate Processor\">" << endl;
581  fout << " <DataArray type=\"Int32\" Name=\"Node\" format=\"ascii\">" << endl;
582  indent = " ";
583  fout << indent;
584  bool localIsGlobal = GlobalOrdinal(nodeMap_->getGlobalNumElements()) == GlobalOrdinal(nodeMap_->getLocalNumElements());
585  for (size_t i = 0; i < uniqueFine.size(); i++) {
586  if (localIsGlobal) {
587  fout << uniqueFine[i] << " "; // if all nodes are on this processor, do not map from local to global
588  } else
589  fout << nodeMap_->getGlobalElement(uniqueFine[i]) << " ";
590  if (i % 10 == 9)
591  fout << endl
592  << indent;
593  }
594  fout << endl;
595  fout << " </DataArray>" << endl;
596  fout << " <DataArray type=\"Int32\" Name=\"Aggregate\" format=\"ascii\">" << endl;
597  fout << indent;
598  for (size_t i = 0; i < uniqueFine.size(); i++) {
599  if (vertex2AggId_[uniqueFine[i]] == -1)
600  fout << vertex2AggId_[uniqueFine[i]] << " ";
601  else
602  fout << aggsOffset_ + vertex2AggId_[uniqueFine[i]] << " ";
603  if (i % 10 == 9)
604  fout << endl
605  << indent;
606  }
607  fout << endl;
608  fout << " </DataArray>" << endl;
609  fout << " <DataArray type=\"Int32\" Name=\"Processor\" format=\"ascii\">" << endl;
610  fout << indent;
611  for (size_t i = 0; i < uniqueFine.size(); i++) {
612  fout << myRank_ << " ";
613  if (i % 20 == 19)
614  fout << endl
615  << indent;
616  }
617  fout << endl;
618  fout << " </DataArray>" << endl;
619  fout << " </PointData>" << endl;
620  fout << " <Points>" << endl;
621  fout << " <DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"ascii\">" << endl;
622  fout << indent;
623  for (size_t i = 0; i < uniqueFine.size(); i++) {
624  fout << xCoords[uniqueFine[i]] << " " << yCoords[uniqueFine[i]] << " ";
625  if (dims_ == 2)
626  fout << "0 ";
627  else
628  fout << zCoords[uniqueFine[i]] << " ";
629  if (i % 3 == 2)
630  fout << endl
631  << indent;
632  }
633  fout << endl;
634  fout << " </DataArray>" << endl;
635  fout << " </Points>" << endl;
636  fout << " <Cells>" << endl;
637  fout << " <DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">" << endl;
638  fout << indent;
639  for (size_t i = 0; i < vertices.size(); i++) {
640  fout << vertices[i] << " ";
641  if (i % 10 == 9)
642  fout << endl
643  << indent;
644  }
645  fout << endl;
646  fout << " </DataArray>" << endl;
647  fout << " <DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">" << endl;
648  fout << indent;
649  int accum = 0;
650  for (size_t i = 0; i < geomSizes.size(); i++) {
651  accum += geomSizes[i];
652  fout << accum << " ";
653  if (i % 10 == 9)
654  fout << endl
655  << indent;
656  }
657  fout << endl;
658  fout << " </DataArray>" << endl;
659  fout << " <DataArray type=\"Int32\" Name=\"types\" format=\"ascii\">" << endl;
660  fout << indent;
661  for (size_t i = 0; i < geomSizes.size(); i++) {
662  switch (geomSizes[i]) {
663  case 1:
664  fout << "1 "; // Point
665  break;
666  case 2:
667  fout << "3 "; // Line
668  break;
669  case 3:
670  fout << "5 "; // Triangle
671  break;
672  default:
673  fout << "7 "; // Polygon
674  }
675  if (i % 30 == 29)
676  fout << endl
677  << indent;
678  }
679  fout << endl;
680  fout << " </DataArray>" << endl;
681  fout << " </Cells>" << endl;
682  fout << " </Piece>" << endl;
683  fout << " </UnstructuredGrid>" << endl;
684  fout << "</VTKFile>" << endl;
685 }
686 
687 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
688 void AggregationExportFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::buildColormap_() const {
689  using namespace std;
690  try {
691  ofstream color("random-colormap.xml");
692  color << "<ColorMap name=\"MueLu-Random\" space=\"RGB\">" << endl;
693  // Give -1, -2, -3 distinctive colors (so that the style functions can have constrasted geometry types)
694  // Do red, orange, yellow to constrast with cool color spectrum for other types
695  color << " <Point x=\"" << CONTRAST_1_ << "\" o=\"1\" r=\"1\" g=\"0\" b=\"0\"/>" << endl;
696  color << " <Point x=\"" << CONTRAST_2_ << "\" o=\"1\" r=\"1\" g=\"0.6\" b=\"0\"/>" << endl;
697  color << " <Point x=\"" << CONTRAST_3_ << "\" o=\"1\" r=\"1\" g=\"1\" b=\"0\"/>" << endl;
698  srand(time(NULL));
699  for (int i = 0; i < 5000; i += 4) {
700  color << " <Point x=\"" << i << "\" o=\"1\" r=\"" << (rand() % 50) / 256.0 << "\" g=\"" << (rand() % 256) / 256.0 << "\" b=\"" << (rand() % 256) / 256.0 << "\"/>" << endl;
701  }
702  color << "</ColorMap>" << endl;
703  color.close();
704  } catch (std::exception& e) {
705  GetOStream(Warnings0) << " Error while building colormap file: " << e.what() << endl;
706  }
707 }
708 
709 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
710 void AggregationExportFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writePVTU_(std::ofstream& pvtu, std::string baseFname, int numProcs) const {
711  using namespace std;
712  // If using vtk, filenameToWrite now contains final, correct ***.vtu filename (for the current proc)
713  // So the root proc will need to use its own filenameToWrite to make a list of the filenames of all other procs to put in
714  // pvtu file.
715  pvtu << "<VTKFile type=\"PUnstructuredGrid\" byte_order=\"LittleEndian\">" << endl;
716  pvtu << " <PUnstructuredGrid GhostLevel=\"0\">" << endl;
717  pvtu << " <PPointData Scalars=\"Node Aggregate Processor\">" << endl;
718  pvtu << " <PDataArray type=\"Int32\" Name=\"Node\"/>" << endl;
719  pvtu << " <PDataArray type=\"Int32\" Name=\"Aggregate\"/>" << endl;
720  pvtu << " <PDataArray type=\"Int32\" Name=\"Processor\"/>" << endl;
721  pvtu << " </PPointData>" << endl;
722  pvtu << " <PPoints>" << endl;
723  pvtu << " <PDataArray type=\"Float64\" NumberOfComponents=\"3\"/>" << endl;
724  pvtu << " </PPoints>" << endl;
725  for (int i = 0; i < numProcs; i++) {
726  // specify the piece for each proc (the replaceAll expression matches what the filenames will be on other procs)
727  pvtu << " <Piece Source=\"" << this->replaceAll(baseFname, "%PROCID", toString(i)) << "\"/>" << endl;
728  }
729  // reference files with graph pieces, if applicable
730  if (doFineGraphEdges_) {
731  for (int i = 0; i < numProcs; i++) {
732  string fn = this->replaceAll(baseFname, "%PROCID", toString(i));
733  pvtu << " <Piece Source=\"" << fn.insert(fn.rfind(".vtu"), "-finegraph") << "\"/>" << endl;
734  }
735  }
736  if (doCoarseGraphEdges_) {
737  for (int i = 0; i < numProcs; i++) {
738  string fn = this->replaceAll(baseFname, "%PROCID", toString(i));
739  pvtu << " <Piece Source=\"" << fn.insert(fn.rfind(".vtu"), "-coarsegraph") << "\"/>" << endl;
740  }
741  }
742  pvtu << " </PUnstructuredGrid>" << endl;
743  pvtu << "</VTKFile>" << endl;
744  pvtu.close();
745 }
746 } // namespace MueLu
747 
748 #endif /* MUELU_AGGREGATIONEXPORTFACTORY_DEF_HPP_ */
MueLu::Warnings0
Important warning messages (one line)
Definition: MueLu_VerbosityLevel.hpp:20
LocalOrdinal
MueLu::DefaultLocalOrdinal LocalOrdinal
Definition: MueLu_UseDefaultTypes.hpp:13
MueLu::toString
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
Definition: MueLu_Utilities_decl.hpp:1102
Xpetra::MultiVectorFactory::Build
static Teuchos::RCP< MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Build(const Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node >> &map, size_t NumVectors, bool zeroOut=true)
GO
GlobalOrdinal GO
Teuchos::ParameterList::get
T & get(const std::string &name, T def_value)
MueLu::FactoryMonitor
Timer to be used in factories. Similar to Monitor but with additional timers.
Definition: MueLu_Monitor.hpp:165
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
MueLu::Runtime0
One-liner description of what is happening.
Definition: MueLu_VerbosityLevel.hpp:25
Teuchos::RCP::get
T * get() const
MueLu_Level.hpp
MueLu::AggregationExportFactory::doConvexHulls
void doConvexHulls(std::vector< int > &vertices, std::vector< int > &geomSizes) const
Definition: MueLu_AggregationExportFactory_def.hpp:303
Teuchos::ParameterList::set
ParameterList & set(std::string const &name, T &&value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Xpetra_Matrix.hpp
MueLu::AggregationExportFactory::writeFile_
void writeFile_(std::ofstream &fout, std::string styleName, std::vector< int > &vertices, std::vector< int > &geomSizes) const
Definition: MueLu_AggregationExportFactory_def.hpp:565
sort
void sort(View &view, const size_t &size)
Teuchos::ParameterList::isParameter
bool isParameter(const std::string &name) const
MueLu::AggregationExportFactory::~AggregationExportFactory
virtual ~AggregationExportFactory()
Destructor.
rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
GlobalOrdinal
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Definition: MueLu_UseDefaultTypes.hpp:14
MueLu::Level
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:63
MueLu::AggregationExportFactory::GetValidParameterList
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Definition: MueLu_AggregationExportFactory_def.hpp:55
MueLu::AggregationExportFactory::writePVTU_
void writePVTU_(std::ofstream &pvtu, std::string baseFname, int numProcs) const
Definition: MueLu_AggregationExportFactory_def.hpp:710
MueLu_AggregationExportFactory_decl.hpp
iter
TransListIter iter
MueLu::AggregationExportFactory::DeclareInput
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
Definition: MueLu_AggregationExportFactory_def.hpp:86
MueLu_Monitor.hpp
MueLu::AggregationExportFactory::doGraphEdges_
void doGraphEdges_(std::ofstream &fout, Teuchos::RCP< Matrix > &A, Teuchos::RCP< LWGraph > &G, bool fine, int dofs) const
Definition: MueLu_AggregationExportFactory_def.hpp:317
MueLu::Level::GetLevelID
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:51
MueLu::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition: MueLu_Exceptions.hpp:36
MueLu::AggregationExportFactory::doJacksPlus_
void doJacksPlus_(std::vector< int > &vertices, std::vector< int > &geomSizes) const
Definition: MueLu_AggregationExportFactory_def.hpp:298
Teuchos::ArrayRCP::get
T * get() const
MueLu::AggregationExportFactory::Build
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
Definition: MueLu_AggregationExportFactory_def.hpp:106
MueLu::replaceAll
void replaceAll(std::string &str, const std::string &from, const std::string &to)
Definition: MueLu_ParameterListUtils.cpp:134
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, const Teuchos::RCP< Teuchos::ParameterList > &plist=Teuchos::null)
Xpetra_ImportFactory.hpp