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