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