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