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 #include "MueLu_Graph.hpp"
64 #include "MueLu_AmalgamationFactory.hpp"
65 #include "MueLu_AmalgamationInfo.hpp"
66 #include "MueLu_Monitor.hpp"
67 #include "MueLu_Utilities.hpp"
68 #include <vector>
69 #include <list>
70 #include <algorithm>
71 #include <string>
72 #include <stdexcept>
73 #include <cstdio>
74 #include <cmath>
75 //For alpha hulls (is optional feature requiring a third-party library)
76 #ifdef HAVE_MUELU_CGAL //Include all headers needed for both 2D and 3D fixed-alpha alpha shapes
77 #include "CGAL/Exact_predicates_inexact_constructions_kernel.h"
78 #include "CGAL/Delaunay_triangulation_2.h"
79 #include "CGAL/Delaunay_triangulation_3.h"
80 #include "CGAL/Alpha_shape_2.h"
81 #include "CGAL/Fixed_alpha_shape_3.h"
82 #include "CGAL/algorithm.h"
83 #endif
84 
85 namespace MueLu {
86 
87  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
88  RCP<const ParameterList> AggregationExportFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetValidParameterList() const {
89  RCP<ParameterList> validParamList = rcp(new ParameterList());
90 
91  std::string output_msg = "Output filename template (%TIMESTEP is replaced by \'Output file: time step\' variable,"
92  "%ITER is replaced by \'Output file: iter\' variable, %LEVELID is replaced level id, %PROCID is replaced by processor id)";
93  std::string output_def = "aggs_level%LEVELID_proc%PROCID.out";
94 
95  validParamList->set< RCP<const FactoryBase> >("A", Teuchos::null, "Factory for A.");
96  validParamList->set< RCP<const FactoryBase> >("Coordinates", Teuchos::null, "Factory for Coordinates.");
97  validParamList->set< RCP<const FactoryBase> >("Graph", Teuchos::null, "Factory for Graph.");
98  validParamList->set< RCP<const FactoryBase> >("Aggregates", Teuchos::null, "Factory for Aggregates.");
99  validParamList->set< RCP<const FactoryBase> >("UnAmalgamationInfo", Teuchos::null, "Factory for UnAmalgamationInfo.");
100  validParamList->set< RCP<const FactoryBase> >("DofsPerNode", Teuchos::null, "Factory for DofsPerNode.");
101  // CMS/BMK: Old style factory-only options. Deprecate me.
102  validParamList->set< std::string > ("Output filename", output_def, output_msg);
103  validParamList->set< int > ("Output file: time step", 0, "time step variable for output file name");
104  validParamList->set< int > ("Output file: iter", 0, "nonlinear iteration variable for output file name");
105 
106  // New-style master list options (here are same defaults as in masterList.xml)
107  validParamList->set< std::string > ("aggregation: output filename", "", "filename for VTK-style visualization output");
108  validParamList->set< int > ("aggregation: output file: time step", 0, "time step variable for output file name");// Remove me?
109  validParamList->set< int > ("aggregation: output file: iter", 0, "nonlinear iteration variable for output file name");//Remove me?
110  validParamList->set<std::string> ("aggregation: output file: agg style", "Point Cloud", "style of aggregate visualization for VTK output");
111  validParamList->set<bool> ("aggregation: output file: fine graph edges", false, "Whether to draw all fine node connections along with the aggregates.");
112  validParamList->set<bool> ("aggregation: output file: coarse graph edges", false, "Whether to draw all coarse node connections along with the aggregates.");
113  validParamList->set<bool> ("aggregation: output file: build colormap", false, "Whether to output a random colormap for ParaView in a separate XML file.");
114  return validParamList;
115  }
116 
117  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
118  void AggregationExportFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::DeclareInput(Level &fineLevel, Level &coarseLevel) const {
119  Input(fineLevel, "Aggregates"); //< factory which created aggregates
120  Input(fineLevel, "DofsPerNode"); //< CoalesceAndDropFactory (needed for DofsPerNode variable)
121  Input(fineLevel, "UnAmalgamationInfo"); //< AmalgamationFactory (needed for UnAmalgamationInfo variable)
122 
123  const ParameterList & pL = GetParameterList();
124  //Only pull in coordinates if the user explicitly requests direct VTK output, so as not to break uses of old code
125  if(pL.isParameter("aggregation: output filename") && pL.get<std::string>("aggregation: output filename").length())
126  {
127  Input(fineLevel, "Coordinates");
128  Input(fineLevel, "A");
129  Input(fineLevel, "Graph");
130  if(pL.get<bool>("aggregation: output file: coarse graph edges"))
131  {
132  Input(coarseLevel, "Coordinates");
133  Input(coarseLevel, "A");
134  Input(coarseLevel, "Graph");
135  }
136  }
137  }
138 
139  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
140  void AggregationExportFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(Level &fineLevel, Level &coarseLevel) const {
141  using namespace std;
142  //Decide which build function to follow, based on input params
143  const ParameterList& pL = GetParameterList();
144  FactoryMonitor m(*this, "AggregationExportFactory", coarseLevel);
145  Teuchos::RCP<Aggregates> aggregates = Get< Teuchos::RCP<Aggregates> >(fineLevel,"Aggregates");
146  Teuchos::RCP<const Teuchos::Comm<int> > comm = aggregates->GetMap()->getComm();
147  int numProcs = comm->getSize();
148  int myRank = comm->getRank();
149  string masterFilename = pL.get<std::string>("aggregation: output filename"); //filename parameter from master list
150  string pvtuFilename = ""; //only root processor will set this
151  string localFilename = pL.get<std::string>("Output filename");
152  string filenameToWrite;
153  bool useVTK = false;
154  doCoarseGraphEdges_ = pL.get<bool>("aggregation: output file: coarse graph edges");
155  doFineGraphEdges_ = pL.get<bool>("aggregation: output file: fine graph edges");
156  if(masterFilename.length())
157  {
158  useVTK = true;
159  filenameToWrite = masterFilename;
160  if(filenameToWrite.rfind(".vtu") == string::npos) //Must have the file extension in the name
161  filenameToWrite.append(".vtu");
162  if(numProcs > 1 && filenameToWrite.rfind("%PROCID") == string::npos) //filename can't be identical between processsors in parallel problem
163  filenameToWrite.insert(filenameToWrite.rfind(".vtu"), "-proc%PROCID");
164  }
165  else
166  filenameToWrite = localFilename;
167  LocalOrdinal DofsPerNode = Get< LocalOrdinal > (fineLevel,"DofsPerNode");
168  Teuchos::RCP<AmalgamationInfo> amalgInfo = Get< RCP<AmalgamationInfo> > (fineLevel,"UnAmalgamationInfo");
169  Teuchos::RCP<Matrix> Amat = Get<RCP<Matrix> >(fineLevel, "A");
170  Teuchos::RCP<Matrix> Ac;
171  if(doCoarseGraphEdges_)
172  Ac = Get<RCP<Matrix> >(coarseLevel, "A");
173  Teuchos::RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType, LocalOrdinal, GlobalOrdinal, Node> > coords = Teuchos::null;
174  Teuchos::RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType, LocalOrdinal, GlobalOrdinal, Node> > coordsCoarse = Teuchos::null;
175  Teuchos::RCP<GraphBase> fineGraph = Teuchos::null;
176  Teuchos::RCP<GraphBase> coarseGraph = Teuchos::null;
177  if(doFineGraphEdges_)
178  fineGraph = Get<RCP<GraphBase> >(fineLevel, "Graph");
179  if(doCoarseGraphEdges_)
180  coarseGraph = Get<RCP<GraphBase> >(coarseLevel, "Graph");
181  if(useVTK) //otherwise leave null, will not be accessed by non-vtk code
182  {
183  coords = Get<RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType, LocalOrdinal, GlobalOrdinal, Node> > >(fineLevel, "Coordinates");
184  if(doCoarseGraphEdges_)
185  coordsCoarse = Get<RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType, LocalOrdinal, GlobalOrdinal, Node> > >(coarseLevel, "Coordinates");
186  dims_ = coords->getNumVectors(); //2D or 3D?
187  if(numProcs > 1)
188  {
189  {
190  RCP<Import> coordImporter = Xpetra::ImportFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(coords->getMap(), Amat->getColMap());
191  RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType, LocalOrdinal, GlobalOrdinal, Node> > ghostedCoords = Xpetra::MultiVectorFactory<typename Teuchos::ScalarTraits<Scalar>::coordinateType, LocalOrdinal, GlobalOrdinal, Node>::Build(Amat->getColMap(), dims_);
192  ghostedCoords->doImport(*coords, *coordImporter, Xpetra::INSERT);
193  coords = ghostedCoords;
194  }
195  if(doCoarseGraphEdges_)
196  {
197  RCP<Import> coordImporter = Xpetra::ImportFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(coordsCoarse->getMap(), Ac->getColMap());
198  RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType, LocalOrdinal, GlobalOrdinal, Node> > ghostedCoords = Xpetra::MultiVectorFactory<typename Teuchos::ScalarTraits<Scalar>::coordinateType, LocalOrdinal, GlobalOrdinal, Node>::Build(Ac->getColMap(), dims_);
199  ghostedCoords->doImport(*coordsCoarse, *coordImporter, Xpetra::INSERT);
200  coordsCoarse = ghostedCoords;
201  }
202  }
203  }
204  GetOStream(Runtime0) << "AggregationExportFactory: DofsPerNode: " << DofsPerNode << std::endl;
205  Teuchos::RCP<LocalOrdinalMultiVector> vertex2AggId_vector = aggregates->GetVertex2AggId();
206  Teuchos::RCP<LocalOrdinalVector> procWinner_vector = aggregates->GetProcWinner();
207  Teuchos::ArrayRCP<LocalOrdinal> vertex2AggId = aggregates->GetVertex2AggId()->getDataNonConst(0);
208  Teuchos::ArrayRCP<LocalOrdinal> procWinner = aggregates->GetProcWinner()->getDataNonConst(0);
209 
210  vertex2AggId_ = vertex2AggId;
211 
212  // prepare for calculating global aggregate ids
213  std::vector<GlobalOrdinal> numAggsGlobal (numProcs, 0);
214  std::vector<GlobalOrdinal> numAggsLocal (numProcs, 0);
215  std::vector<GlobalOrdinal> minGlobalAggId(numProcs, 0);
216 
217  numAggsLocal[myRank] = aggregates->GetNumAggregates();
218  Teuchos::reduceAll(*comm, Teuchos::REDUCE_SUM, numProcs, &numAggsLocal[0], &numAggsGlobal[0]);
219  for (int i = 1; i < Teuchos::as<int>(numAggsGlobal.size()); ++i)
220  {
221  numAggsGlobal [i] += numAggsGlobal[i-1];
222  minGlobalAggId[i] = numAggsGlobal[i-1];
223  }
224  if(numProcs == 0)
225  aggsOffset_ = 0;
226  else
227  aggsOffset_ = minGlobalAggId[myRank];
228  ArrayRCP<LO> aggStart;
229  ArrayRCP<GlobalOrdinal> aggToRowMap;
230  amalgInfo->UnamalgamateAggregates(*aggregates, aggStart, aggToRowMap);
231  int timeStep = pL.get< int > ("Output file: time step");
232  int iter = pL.get< int > ("Output file: iter");
233  filenameToWrite = this->replaceAll(filenameToWrite, "%LEVELID", toString(fineLevel.GetLevelID()));
234  filenameToWrite = this->replaceAll(filenameToWrite, "%TIMESTEP", toString(timeStep));
235  filenameToWrite = this->replaceAll(filenameToWrite, "%ITER", toString(iter));
236  //Proc id MUST be included in vtu filenames to distinguish them (if multiple procs)
237  //In all other cases (else), including processor # in filename is optional
238  string masterStem = "";
239  if(useVTK)
240  {
241  masterStem = filenameToWrite.substr(0, filenameToWrite.rfind(".vtu"));
242  masterStem = this->replaceAll(masterStem, "%PROCID", "");
243  }
244  pvtuFilename = masterStem + "-master.pvtu";
245  string baseFname = filenameToWrite; //get a version of the filename string with the %PROCID token, but without substituting myRank (needed for pvtu output)
246  filenameToWrite = this->replaceAll(filenameToWrite, "%PROCID", toString(myRank));
247  GetOStream(Runtime0) << "AggregationExportFactory: outputfile \"" << filenameToWrite << "\"" << std::endl;
248  ofstream fout(filenameToWrite.c_str());
249  GO numAggs = aggregates->GetNumAggregates();
250  if(!useVTK)
251  {
252  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)
253  GO offset = amalgInfo->GlobalOffset(); // extract offset for global dof ids
254  vector<GlobalOrdinal> nodeIds;
255  for (int i = 0; i < numAggs; ++i) {
256  fout << "Agg " << minGlobalAggId[myRank] + i << " Proc " << myRank << ":";
257 
258  // TODO: Use k+=DofsPerNode instead of ++k and get rid of std::unique call afterwards
259  for (int k = aggStart[i]; k < aggStart[i+1]; ++k) {
260  nodeIds.push_back((aggToRowMap[k] - offset - indexBase) / DofsPerNode + indexBase);
261  }
262 
263  // remove duplicate entries from nodeids
264  std::sort(nodeIds.begin(), nodeIds.end());
265  typename std::vector<GlobalOrdinal>::iterator endLocation = std::unique(nodeIds.begin(), nodeIds.end());
266  nodeIds.erase(endLocation, nodeIds.end());
267 
268  // print out nodeids
269  for(typename std::vector<GlobalOrdinal>::iterator printIt = nodeIds.begin(); printIt != nodeIds.end(); printIt++)
270  fout << " " << *printIt;
271  nodeIds.clear();
272  fout << std::endl;
273  }
274  }
275  else
276  {
277  using namespace std;
278  //Make sure we have coordinates
279  TEUCHOS_TEST_FOR_EXCEPTION(coords.is_null(), Exceptions::RuntimeError,"AggExportFactory could not get coordinates, but they are required for VTK output.");
280  numAggs_ = numAggs;
281  numNodes_ = coords->getLocalLength();
282  //get access to the coord data
283  xCoords_ = Teuchos::arcp_reinterpret_cast<const typename Teuchos::ScalarTraits<Scalar>::coordinateType>(coords->getData(0));
284  yCoords_ = Teuchos::arcp_reinterpret_cast<const typename Teuchos::ScalarTraits<Scalar>::coordinateType>(coords->getData(1));
285  zCoords_ = Teuchos::null;
286  if(doCoarseGraphEdges_)
287  {
288  cx_ = Teuchos::arcp_reinterpret_cast<const typename Teuchos::ScalarTraits<Scalar>::coordinateType>(coordsCoarse->getData(0));
289  cy_ = Teuchos::arcp_reinterpret_cast<const typename Teuchos::ScalarTraits<Scalar>::coordinateType>(coordsCoarse->getData(1));
290  cz_ = Teuchos::null;
291  }
292  if(dims_ == 3)
293  {
294  zCoords_ = Teuchos::arcp_reinterpret_cast<const typename Teuchos::ScalarTraits<Scalar>::coordinateType>(coords->getData(2));
295  if(doCoarseGraphEdges_)
296  cz_ = Teuchos::arcp_reinterpret_cast<const typename Teuchos::ScalarTraits<Scalar>::coordinateType>(coordsCoarse->getData(2));
297  }
298  //Get the sizes of the aggregates to speed up grabbing node IDs
299  aggSizes_ = aggregates->ComputeAggregateSizes();
300  myRank_ = myRank;
301  string aggStyle = "Point Cloud";
302  try
303  {
304  aggStyle = pL.get<string>("aggregation: output file: agg style"); //Let "Point Cloud" be the default style
305  }
306  catch(std::exception& e) {}
307  vector<int> vertices;
308  vector<int> geomSizes;
309  string indent = "";
310  nodeMap_ = Amat->getMap();
311  for(LocalOrdinal i = 0; i < numNodes_; i++)
312  {
313  isRoot_.push_back(aggregates->IsRoot(i));
314  }
315  //If problem is serial and not outputting fine nor coarse graph edges, don't make pvtu file
316  //Otherwise create it
317  if(myRank == 0 && (numProcs != 1 || doCoarseGraphEdges_ || doFineGraphEdges_))
318  {
319  ofstream pvtu(pvtuFilename.c_str());
320  writePVTU_(pvtu, baseFname, numProcs);
321  pvtu.close();
322  }
323  if(aggStyle == "Point Cloud")
324  this->doPointCloud(vertices, geomSizes, numAggs_, numNodes_);
325  else if(aggStyle == "Jacks")
326  this->doJacks(vertices, geomSizes, numAggs_, numNodes_, isRoot_, vertex2AggId_);
327  else if(aggStyle == "Jacks++") //Not actually implemented
328  doJacksPlus_(vertices, geomSizes);
329  else if(aggStyle == "Convex Hulls")
330  doConvexHulls(vertices, geomSizes);
331  else if(aggStyle == "Alpha Hulls")
332  {
333  #ifdef HAVE_MUELU_CGAL
334  doAlphaHulls_(vertices, geomSizes);
335  #else
336  GetOStream(Warnings0) << " Trilinos was not configured with CGAL so Alpha Hulls not available.\n Using Convex Hulls instead." << std::endl;
337  doConvexHulls(vertices, geomSizes);
338  #endif
339  }
340  else
341  {
342  GetOStream(Warnings0) << " Unrecognized agg style.\n Possible values are Point Cloud, Jacks, Jacks++, Convex Hulls and Alpha Hulls.\n Defaulting to Point Cloud." << std::endl;
343  aggStyle = "Point Cloud";
344  this->doPointCloud(vertices, geomSizes, numAggs_, numNodes_);
345  }
346  writeFile_(fout, aggStyle, vertices, geomSizes);
347  if(doCoarseGraphEdges_)
348  {
349  string fname = filenameToWrite;
350  string cEdgeFile = fname.insert(fname.rfind(".vtu"), "-coarsegraph");
351  std::ofstream edgeStream(cEdgeFile.c_str());
352  doGraphEdges_(edgeStream, Ac, coarseGraph, false, DofsPerNode);
353  }
354  if(doFineGraphEdges_)
355  {
356  string fname = filenameToWrite;
357  string fEdgeFile = fname.insert(fname.rfind(".vtu"), "-finegraph");
358  std::ofstream edgeStream(fEdgeFile.c_str());
359  doGraphEdges_(edgeStream, Amat, fineGraph, true, DofsPerNode);
360  }
361  if(myRank == 0 && pL.get<bool>("aggregation: output file: build colormap"))
362  {
363  buildColormap_();
364  }
365  }
366  fout.close();
367  }
368 
369  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
370  void AggregationExportFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::doJacksPlus_(std::vector<int>& /* vertices */, std::vector<int>& /* geomSizes */) const
371  {
372  //TODO
373  }
374 
375  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
376  void AggregationExportFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::doConvexHulls(std::vector<int>& vertices, std::vector<int>& geomSizes) const
377  {
378  if(dims_ == 2)
379  this->doConvexHulls2D(vertices, geomSizes, numAggs_, numNodes_, isRoot_, vertex2AggId_, xCoords_, yCoords_);
380  else
381  this->doConvexHulls3D(vertices, geomSizes, numAggs_, numNodes_, isRoot_, vertex2AggId_, xCoords_, yCoords_, zCoords_);
382  }
383 
384 #ifdef HAVE_MUELU_CGAL
385  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
386  void AggregationExportFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::doAlphaHulls_(std::vector<int>& vertices, std::vector<int>& geomSizes) const
387  {
388  using namespace std;
389  if(dims_ == 2)
390  doAlphaHulls2D_(vertices, geomSizes);
391  else if(dims_ == 3)
392  doAlphaHulls3D_(vertices, geomSizes);
393  }
394 
395  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
396  void AggregationExportFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::doAlphaHulls2D_(std::vector<int>& vertices, std::vector<int>& geomSizes) const
397  {
398  //const double ALPHA_VAL = 2; //Make configurable?
399  using namespace std;
400  //CGAL setup
401  typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
402  typedef K::FT FT;
403  typedef K::Point_2 Point;
404  typedef K::Segment_2 Segment;
405  typedef CGAL::Alpha_shape_vertex_base_2<K> Vb;
406  typedef CGAL::Alpha_shape_face_base_2<K> Fb;
407  typedef CGAL::Triangulation_data_structure_2<Vb,Fb> Tds;
408  typedef CGAL::Delaunay_triangulation_2<K,Tds> Triangulation_2;
409  typedef CGAL::Alpha_shape_2<Triangulation_2> Alpha_shape_2;
410  typedef Alpha_shape_2::Alpha_shape_edges_iterator Alpha_shape_edges_iterator;
411 #if 0 // taw: does not compile with CGAL 4.8
412  for(int i = 0; i < numAggs_; i++)
413  {
414  //Populate a list of Point_2 for this aggregate
415  list<Point> aggPoints;
416  vector<int> aggNodes;
417  for(int j = 0; j < numNodes_; j++)
418  {
419  if(vertex2AggId_[j] == i)
420  {
421  Point p(xCoords_[j], yCoords_[j]);
422  aggPoints.push_back(p);
423  aggNodes.push_back(j);
424  }
425  }
426  Alpha_shape_2 hull(aggPoints.begin(), aggPoints.end(), FT(ALPHA_VAL), Alpha_shape_2::GENERAL);
427  vector<Segment> segments;
428  CGAL::alpha_edges(hull, back_inserter(segments));
429  vertices.reserve(vertices.size() + 2 * segments.size());
430  geomSizes.reserve(geomSizes.size() + segments.size());
431  for(size_t j = 0; j < segments.size(); j++)
432  {
433  for(size_t k = 0; k < aggNodes.size(); k++)
434  {
435  if(fabs(segments[j][0].x == xCoords_[aggNodes[k]]) < 1e-12 && fabs(segments[j][0].y == yCoords_[aggNodes[k]]) < 1e-12)
436  {
437  vertices.push_back(aggNodes[k]);
438  break;
439  }
440  }
441  for(size_t k = 0; k < aggNodes.size(); k++)
442  {
443  if(fabs(segments[j][1].x == xCoords_[aggNodes[k]]) < 1e-12 && fabs(segments[j][1].y == yCoords_[aggNodes[k]]) < 1e-12)
444  {
445  vertices.push_back(aggNodes[k]);
446  break;
447  }
448  }
449  geomSizes.push_back(2); //all cells are line segments
450  }
451  }
452 #endif // if 0
453  }
454 
455  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
456  void AggregationExportFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::doAlphaHulls3D_(std::vector<int>& vertices, std::vector<int>& geomSizes) const
457  {
458  typedef CGAL::Exact_predicates_inexact_constructions_kernel Gt;
459 #if 0 // does not compile with CGAL 4-8
460  typedef CGAL::Alpha_shape_cell_base_3<Gt> Fb;
461  typedef CGAL::Triangulation_data_structure_3<Vb,Fb> Tds;
462  typedef CGAL::Delaunay_triangulation_3<Gt,Tds> Triangulation_3;
463  typedef Gt::Point_3 Point;
464  typedef Alpha_shape_3::Alpha_iterator Alpha_iterator;
465  typedef Alpha_shape_3::Cell_handle Cell_handle;
466  typedef Alpha_shape_3::Vertex_handle Vertex_handle;
467  typedef Alpha_shape_3::Facet Facet;
468  typedef Alpha_shape_3::Edge Edge;
469  typedef Gt::Weighted_point Weighted_point;
470  typedef Gt::Bare_point Bare_point;
471  const double ALPHA_VAL = 2; //Make configurable?
472  using namespace std;
473 
474  for(int i = 0; i < numAggs_; i++)
475  {
476  list<Point> aggPoints;
477  vector<int> aggNodes;
478  for(int j = 0; j < numNodes_; j++)
479  {
480  if(vertex2AggId[j] == i)
481  {
482  Point p(xCoords_[j], yCoords_[j], zCoords_[j]);
483  aggPoints.push_back(p);
484  aggNodes.push_back(j);
485  }
486  }
487  Fixed_alpha_shape_3 hull(aggPoints.begin(), aggPoints.end(), FT(ALPHA_VAL));
488  list<Cell_handle> cells;
489  list<Facet> facets;
490  list<Edge> edges;
491  hull.get_alpha_shape_cells(back_inserter(cells));
492  hull.get_alpha_shape_facets(back_inserter(facets));
493  hull.get_alpha_shape_edges(back_inserter(edges));
494  for(size_t j = 0; j < cells.size(); j++)
495  {
496  Point tetPoints[4];
497  tetPoints[0] = cells[j]->vertex(0);
498  tetPoints[1] = cells[j]->vertex(1);
499  tetPoints[2] = cells[j]->vertex(2);
500  tetPoints[3] = cells[j]->vertex(3);
501  for(int k = 0; k < 4; k++)
502  {
503  for(size_t l = 0; l < aggNodes.size(); l++)
504  {
505  if(fabs(tetPoints[k].x - xCoords_[aggNodes[l]]) < 1e-12 &&
506  fabs(tetPoints[k].y - yCoords_[aggNodes[l]]) < 1e-12 &&
507  fabs(tetPoints[k].z - zCoords_[aggNodes[l]]) < 1e-12)
508  {
509  vertices.push_back(aggNodes[l]);
510  break;
511  }
512  }
513  }
514  geomSizes.push_back(-10); //tetrahedron
515  }
516  for(size_t j = 0; j < facets.size(); j++)
517  {
518  int indices[3];
519  indices[0] = (facets[i].second + 1) % 4;
520  indices[1] = (facets[i].second + 2) % 4;
521  indices[2] = (facets[i].second + 3) % 4;
522  if(facets[i].second % 2 == 0)
523  swap(indices[0], indices[1]);
524  Point facetPts[3];
525  facetPts[0] = facets[i].first->vertex(indices[0])->point();
526  facetPts[1] = facets[i].first->vertex(indices[1])->point();
527  facetPts[2] = facets[i].first->vertex(indices[2])->point();
528  //add triangles in terms of node indices
529  for(size_t l = 0; l < aggNodes.size(); l++)
530  {
531  if(fabs(facetPts[k].x - xCoords_[aggNodes[l]]) < 1e-12 &&
532  fabs(facetPts[k].y - yCoords_[aggNodes[l]]) < 1e-12 &&
533  fabs(facetPts[k].z - zCoords_[aggNodes[l]]) < 1e-12)
534  {
535  vertices.push_back(aggNodes[l]);
536  break;
537  }
538  }
539  geomSizes.push_back(3);
540  }
541  for(size_t j = 0; j < edges.size(); j++)
542  {
543 
544  }
545  }
546 #endif // if 0
547  }
548 #endif
549 
550  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
551  void AggregationExportFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::doGraphEdges_(std::ofstream& fout, Teuchos::RCP<Matrix>& A, Teuchos::RCP<GraphBase>& G, bool fine, int dofs) const
552  {
553  using namespace std;
554  ArrayView<const Scalar> values;
555  ArrayView<const LocalOrdinal> neighbors;
556  //Allow two different colors of connections (by setting "aggregates" scalar to CONTRAST_1 or CONTRAST_2)
557  vector<pair<int, int> > vert1; //vertices (node indices)
558  vector<pair<int, int> > vert2; //size of every cell is assumed to be 2 vertices, since all edges are drawn as lines
559  if(A->isGloballyIndexed())
560  {
561  ArrayView<const GlobalOrdinal> indices;
562  for(GlobalOrdinal globRow = 0; globRow < GlobalOrdinal(A->getGlobalNumRows()); globRow++)
563  {
564  if(dofs == 1)
565  A->getGlobalRowView(globRow, indices, values);
566  neighbors = G->getNeighborVertices((LocalOrdinal) globRow);
567  int gEdge = 0;
568  int aEdge = 0;
569  while(gEdge != int(neighbors.size()))
570  {
571  if(dofs == 1)
572  {
573  if(neighbors[gEdge] == indices[aEdge])
574  {
575  //graph and matrix both have this edge, wasn't filtered, show as color 1
576  vert1.push_back(pair<int, int>(int(globRow), neighbors[gEdge]));
577  gEdge++;
578  aEdge++;
579  }
580  else
581  {
582  //graph contains an edge at gEdge which was filtered from A, show as color 2
583  vert2.push_back(pair<int, int>(int(globRow), neighbors[gEdge]));
584  gEdge++;
585  }
586  }
587  else //for multiple DOF problems, don't try to detect filtered edges and ignore A
588  {
589  vert1.push_back(pair<int, int>(int(globRow), neighbors[gEdge]));
590  gEdge++;
591  }
592  }
593  }
594  }
595  else
596  {
597  ArrayView<const LocalOrdinal> indices;
598  for(LocalOrdinal locRow = 0; locRow < LocalOrdinal(A->getNodeNumRows()); locRow++)
599  {
600  if(dofs == 1)
601  A->getLocalRowView(locRow, indices, values);
602  neighbors = G->getNeighborVertices(locRow);
603  //Add those local indices (columns) to the list of connections (which are pairs of the form (localM, localN))
604  int gEdge = 0;
605  int aEdge = 0;
606  while(gEdge != int(neighbors.size()))
607  {
608  if(dofs == 1)
609  {
610  if(neighbors[gEdge] == indices[aEdge])
611  {
612  vert1.push_back(pair<int, int>(locRow, neighbors[gEdge]));
613  gEdge++;
614  aEdge++;
615  }
616  else
617  {
618  vert2.push_back(pair<int, int>(locRow, neighbors[gEdge]));
619  gEdge++;
620  }
621  }
622  else
623  {
624  vert1.push_back(pair<int, int>(locRow, neighbors[gEdge]));
625  gEdge++;
626  }
627  }
628  }
629  }
630  for(size_t i = 0; i < vert1.size(); i ++)
631  {
632  if(vert1[i].first > vert1[i].second)
633  {
634  int temp = vert1[i].first;
635  vert1[i].first = vert1[i].second;
636  vert1[i].second = temp;
637  }
638  }
639  for(size_t i = 0; i < vert2.size(); i++)
640  {
641  if(vert2[i].first > vert2[i].second)
642  {
643  int temp = vert2[i].first;
644  vert2[i].first = vert2[i].second;
645  vert2[i].second = temp;
646  }
647  }
648  sort(vert1.begin(), vert1.end());
649  vector<pair<int, int> >::iterator newEnd = unique(vert1.begin(), vert1.end()); //remove duplicate edges
650  vert1.erase(newEnd, vert1.end());
651  sort(vert2.begin(), vert2.end());
652  newEnd = unique(vert2.begin(), vert2.end());
653  vert2.erase(newEnd, vert2.end());
654  vector<int> points1;
655  points1.reserve(2 * vert1.size());
656  for(size_t i = 0; i < vert1.size(); i++)
657  {
658  points1.push_back(vert1[i].first);
659  points1.push_back(vert1[i].second);
660  }
661  vector<int> points2;
662  points2.reserve(2 * vert2.size());
663  for(size_t i = 0; i < vert2.size(); i++)
664  {
665  points2.push_back(vert2[i].first);
666  points2.push_back(vert2[i].second);
667  }
668  vector<int> unique1 = this->makeUnique(points1);
669  vector<int> unique2 = this->makeUnique(points2);
670  fout << "<VTKFile type=\"UnstructuredGrid\" byte_order=\"LittleEndian\">" << endl;
671  fout << " <UnstructuredGrid>" << endl;
672  fout << " <Piece NumberOfPoints=\"" << unique1.size() + unique2.size() << "\" NumberOfCells=\"" << vert1.size() + vert2.size() << "\">" << endl;
673  fout << " <PointData Scalars=\"Node Aggregate Processor\">" << endl;
674  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
675  string indent = " ";
676  fout << indent;
677  for(size_t i = 0; i < unique1.size(); i++)
678  {
679  fout << CONTRAST_1_ << " ";
680  if(i % 25 == 24)
681  fout << endl << indent;
682  }
683  for(size_t i = 0; i < unique2.size(); i++)
684  {
685  fout << CONTRAST_2_ << " ";
686  if((i + 2 * vert1.size()) % 25 == 24)
687  fout << endl << indent;
688  }
689  fout << endl;
690  fout << " </DataArray>" << endl;
691  fout << " <DataArray type=\"Int32\" Name=\"Aggregate\" format=\"ascii\">" << endl;
692  fout << indent;
693  for(size_t i = 0; i < unique1.size(); i++)
694  {
695  fout << CONTRAST_1_ << " ";
696  if(i % 25 == 24)
697  fout << endl << indent;
698  }
699  for(size_t i = 0; i < unique2.size(); i++)
700  {
701  fout << CONTRAST_2_ << " ";
702  if((i + 2 * vert2.size()) % 25 == 24)
703  fout << endl << indent;
704  }
705  fout << endl;
706  fout << " </DataArray>" << endl;
707  fout << " <DataArray type=\"Int32\" Name=\"Processor\" format=\"ascii\">" << endl;
708  fout << indent;
709  for(size_t i = 0; i < unique1.size() + unique2.size(); i++)
710  {
711  fout << myRank_ << " ";
712  if(i % 25 == 24)
713  fout << endl << indent;
714  }
715  fout << endl;
716  fout << " </DataArray>" << endl;
717  fout << " </PointData>" << endl;
718  fout << " <Points>" << endl;
719  fout << " <DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"ascii\">" << endl;
720  fout << indent;
721  for(size_t i = 0; i < unique1.size(); i++)
722  {
723  if(fine)
724  {
725  fout << xCoords_[unique1[i]] << " " << yCoords_[unique1[i]] << " ";
726  if(dims_ == 3)
727  fout << zCoords_[unique1[i]] << " ";
728  else
729  fout << "0 ";
730  if(i % 2)
731  fout << endl << indent;
732  }
733  else
734  {
735  fout << cx_[unique1[i]] << " " << cy_[unique1[i]] << " ";
736  if(dims_ == 3)
737  fout << cz_[unique1[i]] << " ";
738  else
739  fout << "0 ";
740  if(i % 2)
741  fout << endl << indent;
742  }
743  }
744  for(size_t i = 0; i < unique2.size(); i++)
745  {
746  if(fine)
747  {
748  fout << xCoords_[unique2[i]] << " " << yCoords_[unique2[i]] << " ";
749  if(dims_ == 3)
750  fout << zCoords_[unique2[i]] << " ";
751  else
752  fout << "0 ";
753  if(i % 2)
754  fout << endl << indent;
755  }
756  else
757  {
758  fout << cx_[unique2[i]] << " " << cy_[unique2[i]] << " ";
759  if(dims_ == 3)
760  fout << cz_[unique2[i]] << " ";
761  else
762  fout << "0 ";
763  if((i + unique1.size()) % 2)
764  fout << endl << indent;
765  }
766  }
767  fout << endl;
768  fout << " </DataArray>" << endl;
769  fout << " </Points>" << endl;
770  fout << " <Cells>" << endl;
771  fout << " <DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">" << endl;
772  fout << indent;
773  for(size_t i = 0; i < points1.size(); i++)
774  {
775  fout << points1[i] << " ";
776  if(i % 10 == 9)
777  fout << endl << indent;
778  }
779  for(size_t i = 0; i < points2.size(); i++)
780  {
781  fout << points2[i] + unique1.size() << " ";
782  if((i + points1.size()) % 10 == 9)
783  fout << endl << indent;
784  }
785  fout << endl;
786  fout << " </DataArray>" << endl;
787  fout << " <DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">" << endl;
788  fout << indent;
789  int offset = 0;
790  for(size_t i = 0; i < vert1.size() + vert2.size(); i++)
791  {
792  offset += 2;
793  fout << offset << " ";
794  if(i % 10 == 9)
795  fout << endl << indent;
796  }
797  fout << endl;
798  fout << " </DataArray>" << endl;
799  fout << " <DataArray type=\"Int32\" Name=\"types\" format=\"ascii\">" << endl;
800  fout << indent;
801  for(size_t i = 0; i < vert1.size() + vert2.size(); i++)
802  {
803  fout << "3 ";
804  if(i % 25 == 24)
805  fout << endl << indent;
806  }
807  fout << endl;
808  fout << " </DataArray>" << endl;
809  fout << " </Cells>" << endl;
810  fout << " </Piece>" << endl;
811  fout << " </UnstructuredGrid>" << endl;
812  fout << "</VTKFile>" << endl;
813  }
814 
815  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
816  void AggregationExportFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeFile_(std::ofstream& fout, std::string styleName, std::vector<int>& vertices, std::vector<int>& geomSizes) const
817  {
818  using namespace std;
819  vector<int> uniqueFine = this->makeUnique(vertices);
820  string indent = " ";
821  fout << "<!--" << styleName << " Aggregates Visualization-->" << endl;
822  fout << "<VTKFile type=\"UnstructuredGrid\" byte_order=\"LittleEndian\">" << endl;
823  fout << " <UnstructuredGrid>" << endl;
824  fout << " <Piece NumberOfPoints=\"" << uniqueFine.size() << "\" NumberOfCells=\"" << geomSizes.size() << "\">" << endl;
825  fout << " <PointData Scalars=\"Node Aggregate Processor\">" << endl;
826  fout << " <DataArray type=\"Int32\" Name=\"Node\" format=\"ascii\">" << endl;
827  indent = " ";
828  fout << indent;
829  bool localIsGlobal = GlobalOrdinal(nodeMap_->getGlobalNumElements()) == GlobalOrdinal(nodeMap_->getNodeNumElements());
830  for(size_t i = 0; i < uniqueFine.size(); i++)
831  {
832  if(localIsGlobal)
833  {
834  fout << uniqueFine[i] << " "; //if all nodes are on this processor, do not map from local to global
835  }
836  else
837  fout << nodeMap_->getGlobalElement(uniqueFine[i]) << " ";
838  if(i % 10 == 9)
839  fout << endl << indent;
840  }
841  fout << endl;
842  fout << " </DataArray>" << endl;
843  fout << " <DataArray type=\"Int32\" Name=\"Aggregate\" format=\"ascii\">" << endl;
844  fout << indent;
845  for(size_t i = 0; i < uniqueFine.size(); i++)
846  {
847  fout << aggsOffset_ + vertex2AggId_[uniqueFine[i]] << " ";
848  if(i % 10 == 9)
849  fout << endl << indent;
850  }
851  fout << endl;
852  fout << " </DataArray>" << endl;
853  fout << " <DataArray type=\"Int32\" Name=\"Processor\" format=\"ascii\">" << endl;
854  fout << indent;
855  for(size_t i = 0; i < uniqueFine.size(); i++)
856  {
857  fout << myRank_ << " ";
858  if(i % 20 == 19)
859  fout << endl << indent;
860  }
861  fout << endl;
862  fout << " </DataArray>" << endl;
863  fout << " </PointData>" << endl;
864  fout << " <Points>" << endl;
865  fout << " <DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"ascii\">" << endl;
866  fout << indent;
867  for(size_t i = 0; i < uniqueFine.size(); i++)
868  {
869  fout << xCoords_[uniqueFine[i]] << " " << yCoords_[uniqueFine[i]] << " ";
870  if(dims_ == 2)
871  fout << "0 ";
872  else
873  fout << zCoords_[uniqueFine[i]] << " ";
874  if(i % 3 == 2)
875  fout << endl << indent;
876  }
877  fout << endl;
878  fout << " </DataArray>" << endl;
879  fout << " </Points>" << endl;
880  fout << " <Cells>" << endl;
881  fout << " <DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">" << endl;
882  fout << indent;
883  for(size_t i = 0; i < vertices.size(); i++)
884  {
885  fout << vertices[i] << " ";
886  if(i % 10 == 9)
887  fout << endl << indent;
888  }
889  fout << endl;
890  fout << " </DataArray>" << endl;
891  fout << " <DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">" << endl;
892  fout << indent;
893  int accum = 0;
894  for(size_t i = 0; i < geomSizes.size(); i++)
895  {
896  accum += geomSizes[i];
897  fout << accum << " ";
898  if(i % 10 == 9)
899  fout << endl << indent;
900  }
901  fout << endl;
902  fout << " </DataArray>" << endl;
903  fout << " <DataArray type=\"Int32\" Name=\"types\" format=\"ascii\">" << endl;
904  fout << indent;
905  for(size_t i = 0; i < geomSizes.size(); i++)
906  {
907  switch(geomSizes[i])
908  {
909  case 1:
910  fout << "1 "; //Point
911  break;
912  case 2:
913  fout << "3 "; //Line
914  break;
915  case 3:
916  fout << "5 "; //Triangle
917  break;
918  default:
919  fout << "7 "; //Polygon
920  }
921  if(i % 30 == 29)
922  fout << endl << indent;
923  }
924  fout << endl;
925  fout << " </DataArray>" << endl;
926  fout << " </Cells>" << endl;
927  fout << " </Piece>" << endl;
928  fout << " </UnstructuredGrid>" << endl;
929  fout << "</VTKFile>" << endl;
930  }
931 
932  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
933  void AggregationExportFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::buildColormap_() const
934  {
935  using namespace std;
936  try
937  {
938  ofstream color("random-colormap.xml");
939  color << "<ColorMap name=\"MueLu-Random\" space=\"RGB\">" << endl;
940  //Give -1, -2, -3 distinctive colors (so that the style functions can have constrasted geometry types)
941  //Do red, orange, yellow to constrast with cool color spectrum for other types
942  color << " <Point x=\"" << CONTRAST_1_ << "\" o=\"1\" r=\"1\" g=\"0\" b=\"0\"/>" << endl;
943  color << " <Point x=\"" << CONTRAST_2_ << "\" o=\"1\" r=\"1\" g=\"0.6\" b=\"0\"/>" << endl;
944  color << " <Point x=\"" << CONTRAST_3_ << "\" o=\"1\" r=\"1\" g=\"1\" b=\"0\"/>" << endl;
945  srand(time(NULL));
946  for(int i = 0; i < 5000; i += 4)
947  {
948  color << " <Point x=\"" << i << "\" o=\"1\" r=\"" << (rand() % 50) / 256.0 << "\" g=\"" << (rand() % 256) / 256.0 << "\" b=\"" << (rand() % 256) / 256.0 << "\"/>" << endl;
949  }
950  color << "</ColorMap>" << endl;
951  color.close();
952  }
953  catch(std::exception& e)
954  {
955  GetOStream(Warnings0) << " Error while building colormap file: " << e.what() << endl;
956  }
957  }
958 
959  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
960  void AggregationExportFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writePVTU_(std::ofstream& pvtu, std::string baseFname, int numProcs) const
961  {
962  using namespace std;
963  //If using vtk, filenameToWrite now contains final, correct ***.vtu filename (for the current proc)
964  //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
965  //pvtu file.
966  pvtu << "<VTKFile type=\"PUnstructuredGrid\" byte_order=\"LittleEndian\">" << endl;
967  pvtu << " <PUnstructuredGrid GhostLevel=\"0\">" << endl;
968  pvtu << " <PPointData Scalars=\"Node Aggregate Processor\">" << endl;
969  pvtu << " <PDataArray type=\"Int32\" Name=\"Node\"/>" << endl;
970  pvtu << " <PDataArray type=\"Int32\" Name=\"Aggregate\"/>" << endl;
971  pvtu << " <PDataArray type=\"Int32\" Name=\"Processor\"/>" << endl;
972  pvtu << " </PPointData>" << endl;
973  pvtu << " <PPoints>" << endl;
974  pvtu << " <PDataArray type=\"Float64\" NumberOfComponents=\"3\"/>" << endl;
975  pvtu << " </PPoints>" << endl;
976  for(int i = 0; i < numProcs; i++)
977  {
978  //specify the piece for each proc (the replaceAll expression matches what the filenames will be on other procs)
979  pvtu << " <Piece Source=\"" << this->replaceAll(baseFname, "%PROCID", toString(i)) << "\"/>" << endl;
980  }
981  //reference files with graph pieces, if applicable
982  if(doFineGraphEdges_)
983  {
984  for(int i = 0; i < numProcs; i++)
985  {
986  string fn = this->replaceAll(baseFname, "%PROCID", toString(i));
987  pvtu << " <Piece Source=\"" << fn.insert(fn.rfind(".vtu"), "-finegraph") << "\"/>" << endl;
988  }
989  }
990  if(doCoarseGraphEdges_)
991  {
992  for(int i = 0; i < numProcs; i++)
993  {
994  string fn = this->replaceAll(baseFname, "%PROCID", toString(i));
995  pvtu << " <Piece Source=\"" << fn.insert(fn.rfind(".vtu"), "-coarsegraph") << "\"/>" << endl;
996  }
997  }
998  pvtu << " </PUnstructuredGrid>" << endl;
999  pvtu << "</VTKFile>" << endl;
1000  pvtu.close();
1001  }
1002 } // namespace MueLu
1003 
1004 #endif /* MUELU_AGGREGATIONEXPORTFACTORY_DEF_HPP_ */
MueLu::Warnings0
Important warning messages (one line)
Definition: MueLu_VerbosityLevel.hpp:57
Xpetra::ImportFactory::Build
static RCP< Import< LocalOrdinal, GlobalOrdinal, Node > > Build(const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &source, const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &target)
MueLu::AggregationExportFactory::doAlphaHulls3D_
void doAlphaHulls3D_(std::vector< int > &vertices, std::vector< int > &geomSizes) const
Definition: MueLu_AggregationExportFactory_def.hpp:456
MueLu::toString
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
Definition: MueLu_Utilities_decl.hpp:990
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:209
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
MueLu::AggregationExportFactory::doAlphaHulls_
void doAlphaHulls_(std::vector< int > &vertices, std::vector< int > &geomSizes) const
Definition: MueLu_AggregationExportFactory_def.hpp:386
MueLu::Runtime0
One-liner description of what is happening.
Definition: MueLu_VerbosityLevel.hpp:62
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:376
Xpetra_Matrix.hpp
MueLu::AggregationExportFactory::doGraphEdges_
void doGraphEdges_(std::ofstream &fout, Teuchos::RCP< Matrix > &A, Teuchos::RCP< GraphBase > &G, bool fine, int dofs) const
Definition: MueLu_AggregationExportFactory_def.hpp:551
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:816
Teuchos::ParameterList::isParameter
bool isParameter(const std::string &name) const
Xpetra_MultiVectorFactory.hpp
rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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:88
MueLu::AggregationExportFactory::writePVTU_
void writePVTU_(std::ofstream &pvtu, std::string baseFname, int numProcs) const
Definition: MueLu_AggregationExportFactory_def.hpp:960
Teuchos::reduceAll
TEUCHOS_DEPRECATED void reduceAll(const Comm< Ordinal > &comm, const EReductionType reductType, const Packet &send, Packet *globalReduct)
MueLu_AggregationExportFactory_decl.hpp
iter
TransListIter iter
MueLu::AggregationExportFactory::DeclareInput
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
Definition: MueLu_AggregationExportFactory_def.hpp:118
MueLu_Monitor.hpp
MueLu::Level::GetLevelID
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:76
Xpetra_CrsMatrixWrap.hpp
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:370
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:140
MueLu::replaceAll
void replaceAll(std::string &str, const std::string &from, const std::string &to)
Definition: MueLu_ParameterListUtils.cpp:175
MueLu::AggregationExportFactory::doAlphaHulls2D_
void doAlphaHulls2D_(std::vector< int > &vertices, std::vector< int > &geomSizes) const
Definition: MueLu_AggregationExportFactory_def.hpp:396
Xpetra_ImportFactory.hpp