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<double, LocalOrdinal, GlobalOrdinal, Node> > coords = Teuchos::null;
174  Teuchos::RCP<Xpetra::MultiVector<double, 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<double, LocalOrdinal, GlobalOrdinal, Node> > >(fineLevel, "Coordinates");
184  if(doCoarseGraphEdges_)
185  coordsCoarse = Get<RCP<Xpetra::MultiVector<double, 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<double, LocalOrdinal, GlobalOrdinal, Node> > ghostedCoords = Xpetra::MultiVectorFactory<double, 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<double, LocalOrdinal, GlobalOrdinal, Node> > ghostedCoords = Xpetra::MultiVectorFactory<double, 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<LocalOrdinalVector> 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 double>(coords->getData(0));
284  yCoords_ = Teuchos::arcp_reinterpret_cast<const double>(coords->getData(1));
285  zCoords_ = Teuchos::null;
286  if(doCoarseGraphEdges_)
287  {
288  cx_ = Teuchos::arcp_reinterpret_cast<const double>(coordsCoarse->getData(0));
289  cy_ = Teuchos::arcp_reinterpret_cast<const double>(coordsCoarse->getData(1));
290  cz_ = Teuchos::null;
291  }
292  if(dims_ == 3)
293  {
294  zCoords_ = Teuchos::arcp_reinterpret_cast<const double>(coords->getData(2));
295  if(doCoarseGraphEdges_)
296  cz_ = Teuchos::arcp_reinterpret_cast<const double>(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(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_, zCoords_);
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  for(int i = 0; i < numAggs_; i++)
412  {
413  //Populate a list of Point_2 for this aggregate
414  list<Point> aggPoints;
415  vector<int> aggNodes;
416  for(int j = 0; j < numNodes_; j++)
417  {
418  if(vertex2AggId_[j] == i)
419  {
420  Point p(xCoords_[j], yCoords_[j]);
421  aggPoints.push_back(p);
422  aggNodes.push_back(j);
423  }
424  }
425  Alpha_shape_2 hull(aggPoints.begin(), aggPoints.end(), FT(ALPHA_VAL), Alpha_shape_2::GENERAL);
426  vector<Segment> segments;
427  alpha_edges(hull, back_inserter(segments));
428  vertices.reserve(vertices.size() + 2 * segments.size());
429  geomSizes.reserve(geomSizes.size() + segments.size());
430  for(size_t j = 0; j < segments.size(); j++)
431  {
432  for(size_t k = 0; k < aggNodes.size(); k++)
433  {
434  if(fabs(segments[j][0].x == xCoords_[aggNodes[k]]) < 1e-12 && fabs(segments[j][0].y == yCoords_[aggNodes[k]]) < 1e-12)
435  {
436  vertices.push_back(aggNodes[k]);
437  break;
438  }
439  }
440  for(size_t k = 0; k < aggNodes.size(); k++)
441  {
442  if(fabs(segments[j][1].x == xCoords_[aggNodes[k]]) < 1e-12 && fabs(segments[j][1].y == yCoords_[aggNodes[k]]) < 1e-12)
443  {
444  vertices.push_back(aggNodes[k]);
445  break;
446  }
447  }
448  geomSizes.push_back(2); //all cells are line segments
449  }
450  }
451  }
452 
453  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
454  void AggregationExportFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::doAlphaHulls3D_(std::vector<int>& vertices, std::vector<int>& geomSizes) const
455  {
456  typedef CGAL::Exact_predicates_inexact_constructions_kernel Gt;
457  typedef CGAL::Alpha_shape_cell_base_3<Gt> Fb;
458  typedef CGAL::Triangulation_data_structure_3<Vb,Fb> Tds;
459  typedef CGAL::Delaunay_triangulation_3<Gt,Tds> Triangulation_3;
460  typedef Gt::Point_3 Point;
461  typedef Alpha_shape_3::Alpha_iterator Alpha_iterator;
462  typedef Alpha_shape_3::Cell_handle Cell_handle;
463  typedef Alpha_shape_3::Vertex_handle Vertex_handle;
464  typedef Alpha_shape_3::Facet Facet;
465  typedef Alpha_shape_3::Edge Edge;
466  typedef Gt::Weighted_point Weighted_point;
467  typedef Gt::Bare_point Bare_point;
468  const double ALPHA_VAL = 2; //Make configurable?
469  using namespace std;
470  for(int i = 0; i < numAggs_; i++)
471  {
472  list<Point> aggPoints;
473  vector<int> aggNodes;
474  for(int j = 0; j < numNodes_; j++)
475  {
476  if(vertex2AggId[j] == i)
477  {
478  Point p(xCoords_[j], yCoords_[j], zCoords_[j]);
479  aggPoints.push_back(p);
480  aggNodes.push_back(j);
481  }
482  }
483  Fixed_alpha_shape_3 hull(aggPoints.begin(), aggPoints.end(), FT(ALPHA_VAL));
484  list<Cell_handle> cells;
485  list<Facet> facets;
486  list<Edge> edges;
487  hull.get_alpha_shape_cells(back_inserter(cells));
488  hull.get_alpha_shape_facets(back_inserter(facets));
489  hull.get_alpha_shape_edges(back_inserter(edges));
490  for(size_t j = 0; j < cells.size(); j++)
491  {
492  Point tetPoints[4];
493  tetPoints[0] = cells[j]->vertex(0);
494  tetPoints[1] = cells[j]->vertex(1);
495  tetPoints[2] = cells[j]->vertex(2);
496  tetPoints[3] = cells[j]->vertex(3);
497  for(int k = 0; k < 4; k++)
498  {
499  for(size_t l = 0; l < aggNodes.size(); l++)
500  {
501  if(fabs(tetPoints[k].x - xCoords_[aggNodes[l]]) < 1e-12 &&
502  fabs(tetPoints[k].y - yCoords_[aggNodes[l]]) < 1e-12 &&
503  fabs(tetPoints[k].z - zCoords_[aggNodes[l]]) < 1e-12)
504  {
505  vertices.push_back(aggNodes[l]);
506  break;
507  }
508  }
509  }
510  geomSizes.push_back(-10); //tetrahedron
511  }
512  for(size_t j = 0; j < facets.size(); j++)
513  {
514  int indices[3];
515  indices[0] = (facets[i].second + 1) % 4;
516  indices[1] = (facets[i].second + 2) % 4;
517  indices[2] = (facets[i].second + 3) % 4;
518  if(facets[i].second % 2 == 0)
519  swap(indices[0], indices[1]);
520  Point facetPts[3];
521  facetPts[0] = facets[i].first->vertex(indices[0])->point();
522  facetPts[1] = facets[i].first->vertex(indices[1])->point();
523  facetPts[2] = facets[i].first->vertex(indices[2])->point();
524  //add triangles in terms of node indices
525  for(size_t l = 0; l < aggNodes.size(); l++)
526  {
527  if(fabs(facetPts[k].x - xCoords_[aggNodes[l]]) < 1e-12 &&
528  fabs(facetPts[k].y - yCoords_[aggNodes[l]]) < 1e-12 &&
529  fabs(facetPts[k].z - zCoords_[aggNodes[l]]) < 1e-12)
530  {
531  vertices.push_back(aggNodes[l]);
532  break;
533  }
534  }
535  geomSizes.push_back(3);
536  }
537  for(size_t j = 0; j < edges.size(); j++)
538  {
539 
540  }
541  }
542  }
543 #endif
544 
545  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
546  void AggregationExportFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::doGraphEdges_(std::ofstream& fout, Teuchos::RCP<Matrix>& A, Teuchos::RCP<GraphBase>& G, bool fine, int dofs) const
547  {
548  using namespace std;
549  ArrayView<const Scalar> values;
550  ArrayView<const LocalOrdinal> neighbors;
551  //Allow two different colors of connections (by setting "aggregates" scalar to CONTRAST_1 or CONTRAST_2)
552  vector<pair<int, int> > vert1; //vertices (node indices)
553  vector<pair<int, int> > vert2; //size of every cell is assumed to be 2 vertices, since all edges are drawn as lines
554  if(A->isGloballyIndexed())
555  {
556  ArrayView<const GlobalOrdinal> indices;
557  for(GlobalOrdinal globRow = 0; globRow < GlobalOrdinal(A->getGlobalNumRows()); globRow++)
558  {
559  if(dofs == 1)
560  A->getGlobalRowView(globRow, indices, values);
561  neighbors = G->getNeighborVertices((LocalOrdinal) globRow);
562  int gEdge = 0;
563  int aEdge = 0;
564  while(gEdge != int(neighbors.size()))
565  {
566  if(dofs == 1)
567  {
568  if(neighbors[gEdge] == indices[aEdge])
569  {
570  //graph and matrix both have this edge, wasn't filtered, show as color 1
571  vert1.push_back(pair<int, int>(int(globRow), neighbors[gEdge]));
572  gEdge++;
573  aEdge++;
574  }
575  else
576  {
577  //graph contains an edge at gEdge which was filtered from A, show as color 2
578  vert2.push_back(pair<int, int>(int(globRow), neighbors[gEdge]));
579  gEdge++;
580  }
581  }
582  else //for multiple DOF problems, don't try to detect filtered edges and ignore A
583  {
584  vert1.push_back(pair<int, int>(int(globRow), neighbors[gEdge]));
585  gEdge++;
586  }
587  }
588  }
589  }
590  else
591  {
592  ArrayView<const LocalOrdinal> indices;
593  for(LocalOrdinal locRow = 0; locRow < LocalOrdinal(A->getNodeNumRows()); locRow++)
594  {
595  if(dofs == 1)
596  A->getLocalRowView(locRow, indices, values);
597  neighbors = G->getNeighborVertices(locRow);
598  //Add those local indices (columns) to the list of connections (which are pairs of the form (localM, localN))
599  int gEdge = 0;
600  int aEdge = 0;
601  while(gEdge != int(neighbors.size()))
602  {
603  if(dofs == 1)
604  {
605  if(neighbors[gEdge] == indices[aEdge])
606  {
607  vert1.push_back(pair<int, int>(locRow, neighbors[gEdge]));
608  gEdge++;
609  aEdge++;
610  }
611  else
612  {
613  vert2.push_back(pair<int, int>(locRow, neighbors[gEdge]));
614  gEdge++;
615  }
616  }
617  else
618  {
619  vert1.push_back(pair<int, int>(locRow, neighbors[gEdge]));
620  gEdge++;
621  }
622  }
623  }
624  }
625  for(size_t i = 0; i < vert1.size(); i ++)
626  {
627  if(vert1[i].first > vert1[i].second)
628  {
629  int temp = vert1[i].first;
630  vert1[i].first = vert1[i].second;
631  vert1[i].second = temp;
632  }
633  }
634  for(size_t i = 0; i < vert2.size(); i++)
635  {
636  if(vert2[i].first > vert2[i].second)
637  {
638  int temp = vert2[i].first;
639  vert2[i].first = vert2[i].second;
640  vert2[i].second = temp;
641  }
642  }
643  sort(vert1.begin(), vert1.end());
644  vector<pair<int, int> >::iterator newEnd = unique(vert1.begin(), vert1.end()); //remove duplicate edges
645  vert1.erase(newEnd, vert1.end());
646  sort(vert2.begin(), vert2.end());
647  newEnd = unique(vert2.begin(), vert2.end());
648  vert2.erase(newEnd, vert2.end());
649  vector<int> points1;
650  points1.reserve(2 * vert1.size());
651  for(size_t i = 0; i < vert1.size(); i++)
652  {
653  points1.push_back(vert1[i].first);
654  points1.push_back(vert1[i].second);
655  }
656  vector<int> points2;
657  points2.reserve(2 * vert2.size());
658  for(size_t i = 0; i < vert2.size(); i++)
659  {
660  points2.push_back(vert2[i].first);
661  points2.push_back(vert2[i].second);
662  }
663  vector<int> unique1 = this->makeUnique(points1);
664  vector<int> unique2 = this->makeUnique(points2);
665  fout << "<VTKFile type=\"UnstructuredGrid\" byte_order=\"LittleEndian\">" << endl;
666  fout << " <UnstructuredGrid>" << endl;
667  fout << " <Piece NumberOfPoints=\"" << unique1.size() + unique2.size() << "\" NumberOfCells=\"" << vert1.size() + vert2.size() << "\">" << endl;
668  fout << " <PointData Scalars=\"Node Aggregate Processor\">" << endl;
669  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
670  string indent = " ";
671  fout << indent;
672  for(size_t i = 0; i < unique1.size(); i++)
673  {
674  fout << CONTRAST_1_ << " ";
675  if(i % 25 == 24)
676  fout << endl << indent;
677  }
678  for(size_t i = 0; i < unique2.size(); i++)
679  {
680  fout << CONTRAST_2_ << " ";
681  if((i + 2 * vert1.size()) % 25 == 24)
682  fout << endl << indent;
683  }
684  fout << endl;
685  fout << " </DataArray>" << endl;
686  fout << " <DataArray type=\"Int32\" Name=\"Aggregate\" format=\"ascii\">" << endl;
687  fout << indent;
688  for(size_t i = 0; i < unique1.size(); i++)
689  {
690  fout << CONTRAST_1_ << " ";
691  if(i % 25 == 24)
692  fout << endl << indent;
693  }
694  for(size_t i = 0; i < unique2.size(); i++)
695  {
696  fout << CONTRAST_2_ << " ";
697  if((i + 2 * vert2.size()) % 25 == 24)
698  fout << endl << indent;
699  }
700  fout << endl;
701  fout << " </DataArray>" << endl;
702  fout << " <DataArray type=\"Int32\" Name=\"Processor\" format=\"ascii\">" << endl;
703  fout << indent;
704  for(size_t i = 0; i < unique1.size() + unique2.size(); i++)
705  {
706  fout << myRank_ << " ";
707  if(i % 25 == 24)
708  fout << endl << indent;
709  }
710  fout << endl;
711  fout << " </DataArray>" << endl;
712  fout << " </PointData>" << endl;
713  fout << " <Points>" << endl;
714  fout << " <DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"ascii\">" << endl;
715  fout << indent;
716  for(size_t i = 0; i < unique1.size(); i++)
717  {
718  if(fine)
719  {
720  fout << xCoords_[unique1[i]] << " " << yCoords_[unique1[i]] << " ";
721  if(dims_ == 3)
722  fout << zCoords_[unique1[i]] << " ";
723  else
724  fout << "0 ";
725  if(i % 2)
726  fout << endl << indent;
727  }
728  else
729  {
730  fout << cx_[unique1[i]] << " " << cy_[unique1[i]] << " ";
731  if(dims_ == 3)
732  fout << cz_[unique1[i]] << " ";
733  else
734  fout << "0 ";
735  if(i % 2)
736  fout << endl << indent;
737  }
738  }
739  for(size_t i = 0; i < unique2.size(); i++)
740  {
741  if(fine)
742  {
743  fout << xCoords_[unique2[i]] << " " << yCoords_[unique2[i]] << " ";
744  if(dims_ == 3)
745  fout << zCoords_[unique2[i]] << " ";
746  else
747  fout << "0 ";
748  if(i % 2)
749  fout << endl << indent;
750  }
751  else
752  {
753  fout << cx_[unique2[i]] << " " << cy_[unique2[i]] << " ";
754  if(dims_ == 3)
755  fout << cz_[unique2[i]] << " ";
756  else
757  fout << "0 ";
758  if((i + unique1.size()) % 2)
759  fout << endl << indent;
760  }
761  }
762  fout << endl;
763  fout << " </DataArray>" << endl;
764  fout << " </Points>" << endl;
765  fout << " <Cells>" << endl;
766  fout << " <DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">" << endl;
767  fout << indent;
768  for(size_t i = 0; i < points1.size(); i++)
769  {
770  fout << points1[i] << " ";
771  if(i % 10 == 9)
772  fout << endl << indent;
773  }
774  for(size_t i = 0; i < points2.size(); i++)
775  {
776  fout << points2[i] + unique1.size() << " ";
777  if((i + points1.size()) % 10 == 9)
778  fout << endl << indent;
779  }
780  fout << endl;
781  fout << " </DataArray>" << endl;
782  fout << " <DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">" << endl;
783  fout << indent;
784  int offset = 0;
785  for(size_t i = 0; i < vert1.size() + vert2.size(); i++)
786  {
787  offset += 2;
788  fout << offset << " ";
789  if(i % 10 == 9)
790  fout << endl << indent;
791  }
792  fout << endl;
793  fout << " </DataArray>" << endl;
794  fout << " <DataArray type=\"Int32\" Name=\"types\" format=\"ascii\">" << endl;
795  fout << indent;
796  for(size_t i = 0; i < vert1.size() + vert2.size(); i++)
797  {
798  fout << "3 ";
799  if(i % 25 == 24)
800  fout << endl << indent;
801  }
802  fout << endl;
803  fout << " </DataArray>" << endl;
804  fout << " </Cells>" << endl;
805  fout << " </Piece>" << endl;
806  fout << " </UnstructuredGrid>" << endl;
807  fout << "</VTKFile>" << endl;
808  }
809 
810  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
811  void AggregationExportFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeFile_(std::ofstream& fout, std::string styleName, std::vector<int>& vertices, std::vector<int>& geomSizes) const
812  {
813  using namespace std;
814  vector<int> uniqueFine = this->makeUnique(vertices);
815  string indent = " ";
816  fout << "<!--" << styleName << " Aggregates Visualization-->" << endl;
817  fout << "<VTKFile type=\"UnstructuredGrid\" byte_order=\"LittleEndian\">" << endl;
818  fout << " <UnstructuredGrid>" << endl;
819  fout << " <Piece NumberOfPoints=\"" << uniqueFine.size() << "\" NumberOfCells=\"" << geomSizes.size() << "\">" << endl;
820  fout << " <PointData Scalars=\"Node Aggregate Processor\">" << endl;
821  fout << " <DataArray type=\"Int32\" Name=\"Node\" format=\"ascii\">" << endl;
822  indent = " ";
823  fout << indent;
824  bool localIsGlobal = GlobalOrdinal(nodeMap_->getGlobalNumElements()) == GlobalOrdinal(nodeMap_->getNodeNumElements());
825  for(size_t i = 0; i < uniqueFine.size(); i++)
826  {
827  if(localIsGlobal)
828  {
829  fout << uniqueFine[i] << " "; //if all nodes are on this processor, do not map from local to global
830  }
831  else
832  fout << nodeMap_->getGlobalElement(uniqueFine[i]) << " ";
833  if(i % 10 == 9)
834  fout << endl << indent;
835  }
836  fout << endl;
837  fout << " </DataArray>" << endl;
838  fout << " <DataArray type=\"Int32\" Name=\"Aggregate\" format=\"ascii\">" << endl;
839  fout << indent;
840  for(size_t i = 0; i < uniqueFine.size(); i++)
841  {
842  fout << aggsOffset_ + vertex2AggId_[uniqueFine[i]] << " ";
843  if(i % 10 == 9)
844  fout << endl << indent;
845  }
846  fout << endl;
847  fout << " </DataArray>" << endl;
848  fout << " <DataArray type=\"Int32\" Name=\"Processor\" format=\"ascii\">" << endl;
849  fout << indent;
850  for(size_t i = 0; i < uniqueFine.size(); i++)
851  {
852  fout << myRank_ << " ";
853  if(i % 20 == 19)
854  fout << endl << indent;
855  }
856  fout << endl;
857  fout << " </DataArray>" << endl;
858  fout << " </PointData>" << endl;
859  fout << " <Points>" << endl;
860  fout << " <DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"ascii\">" << endl;
861  fout << indent;
862  for(size_t i = 0; i < uniqueFine.size(); i++)
863  {
864  fout << xCoords_[uniqueFine[i]] << " " << yCoords_[uniqueFine[i]] << " ";
865  if(dims_ == 2)
866  fout << "0 ";
867  else
868  fout << zCoords_[uniqueFine[i]] << " ";
869  if(i % 3 == 2)
870  fout << endl << indent;
871  }
872  fout << endl;
873  fout << " </DataArray>" << endl;
874  fout << " </Points>" << endl;
875  fout << " <Cells>" << endl;
876  fout << " <DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">" << endl;
877  fout << indent;
878  for(size_t i = 0; i < vertices.size(); i++)
879  {
880  fout << vertices[i] << " ";
881  if(i % 10 == 9)
882  fout << endl << indent;
883  }
884  fout << endl;
885  fout << " </DataArray>" << endl;
886  fout << " <DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">" << endl;
887  fout << indent;
888  int accum = 0;
889  for(size_t i = 0; i < geomSizes.size(); i++)
890  {
891  accum += geomSizes[i];
892  fout << accum << " ";
893  if(i % 10 == 9)
894  fout << endl << indent;
895  }
896  fout << endl;
897  fout << " </DataArray>" << endl;
898  fout << " <DataArray type=\"Int32\" Name=\"types\" format=\"ascii\">" << endl;
899  fout << indent;
900  for(size_t i = 0; i < geomSizes.size(); i++)
901  {
902  switch(geomSizes[i])
903  {
904  case 1:
905  fout << "1 "; //Point
906  break;
907  case 2:
908  fout << "3 "; //Line
909  break;
910  case 3:
911  fout << "5 "; //Triangle
912  break;
913  default:
914  fout << "7 "; //Polygon
915  }
916  if(i % 30 == 29)
917  fout << endl << indent;
918  }
919  fout << endl;
920  fout << " </DataArray>" << endl;
921  fout << " </Cells>" << endl;
922  fout << " </Piece>" << endl;
923  fout << " </UnstructuredGrid>" << endl;
924  fout << "</VTKFile>" << endl;
925  }
926 
927  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
929  {
930  using namespace std;
931  try
932  {
933  ofstream color("random-colormap.xml");
934  color << "<ColorMap name=\"MueLu-Random\" space=\"RGB\">" << endl;
935  //Give -1, -2, -3 distinctive colors (so that the style functions can have constrasted geometry types)
936  //Do red, orange, yellow to constrast with cool color spectrum for other types
937  color << " <Point x=\"" << CONTRAST_1_ << "\" o=\"1\" r=\"1\" g=\"0\" b=\"0\"/>" << endl;
938  color << " <Point x=\"" << CONTRAST_2_ << "\" o=\"1\" r=\"1\" g=\"0.6\" b=\"0\"/>" << endl;
939  color << " <Point x=\"" << CONTRAST_3_ << "\" o=\"1\" r=\"1\" g=\"1\" b=\"0\"/>" << endl;
940  srand(time(NULL));
941  for(int i = 0; i < 5000; i += 4)
942  {
943  color << " <Point x=\"" << i << "\" o=\"1\" r=\"" << (rand() % 50) / 256.0 << "\" g=\"" << (rand() % 256) / 256.0 << "\" b=\"" << (rand() % 256) / 256.0 << "\"/>" << endl;
944  }
945  color << "</ColorMap>" << endl;
946  color.close();
947  }
948  catch(exception& e)
949  {
950  GetOStream(Warnings0) << " Error while building colormap file: " << e.what() << endl;
951  }
952  }
953 
954  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
955  void AggregationExportFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writePVTU_(std::ofstream& pvtu, std::string baseFname, int numProcs) const
956  {
957  using namespace std;
958  //If using vtk, filenameToWrite now contains final, correct ***.vtu filename (for the current proc)
959  //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
960  //pvtu file.
961  pvtu << "<VTKFile type=\"PUnstructuredGrid\" byte_order=\"LittleEndian\">" << endl;
962  pvtu << " <PUnstructuredGrid GhostLevel=\"0\">" << endl;
963  pvtu << " <PPointData Scalars=\"Node Aggregate Processor\">" << endl;
964  pvtu << " <PDataArray type=\"Int32\" Name=\"Node\"/>" << endl;
965  pvtu << " <PDataArray type=\"Int32\" Name=\"Aggregate\"/>" << endl;
966  pvtu << " <PDataArray type=\"Int32\" Name=\"Processor\"/>" << endl;
967  pvtu << " </PPointData>" << endl;
968  pvtu << " <PPoints>" << endl;
969  pvtu << " <PDataArray type=\"Float64\" NumberOfComponents=\"3\"/>" << endl;
970  pvtu << " </PPoints>" << endl;
971  for(int i = 0; i < numProcs; i++)
972  {
973  //specify the piece for each proc (the replaceAll expression matches what the filenames will be on other procs)
974  pvtu << " <Piece Source=\"" << this->replaceAll(baseFname, "%PROCID", toString(i)) << "\"/>" << endl;
975  }
976  //reference files with graph pieces, if applicable
977  if(doFineGraphEdges_)
978  {
979  for(int i = 0; i < numProcs; i++)
980  {
981  string fn = this->replaceAll(baseFname, "%PROCID", toString(i));
982  pvtu << " <Piece Source=\"" << fn.insert(fn.rfind(".vtu"), "-finegraph") << "\"/>" << endl;
983  }
984  }
985  if(doCoarseGraphEdges_)
986  {
987  for(int i = 0; i < numProcs; i++)
988  {
989  string fn = this->replaceAll(baseFname, "%PROCID", toString(i));
990  pvtu << " <Piece Source=\"" << fn.insert(fn.rfind(".vtu"), "-coarsegraph") << "\"/>" << endl;
991  }
992  }
993  pvtu << " </PUnstructuredGrid>" << endl;
994  pvtu << "</VTKFile>" << endl;
995  pvtu.close();
996  }
997 } // namespace MueLu
998 
999 #endif /* MUELU_AGGREGATIONEXPORTFACTORY_DEF_HPP_ */
Important warning messages (one line)
void doAlphaHulls3D_(std::vector< int > &vertices, std::vector< int > &geomSizes) const
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
Timer to be used in factories. Similar to Monitor but with additional timers.
void doAlphaHulls_(std::vector< int > &vertices, std::vector< int > &geomSizes) const
One-liner description of what is happening.
void doConvexHulls(std::vector< int > &vertices, std::vector< int > &geomSizes) const
Namespace for MueLu classes and methods.
void doGraphEdges_(std::ofstream &fout, Teuchos::RCP< Matrix > &A, Teuchos::RCP< GraphBase > &G, bool fine, int dofs) const
void writeFile_(std::ofstream &fout, std::string styleName, std::vector< int > &vertices, std::vector< int > &geomSizes) const
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
void writePVTU_(std::ofstream &pvtu, std::string baseFname, int numProcs) const
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
Exception throws to report errors in the internal logical of the program.
void doJacksPlus_(std::vector< int > &vertices, std::vector< int > &geomSizes) const
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
void replaceAll(std::string &str, const std::string &from, const std::string &to)
void doAlphaHulls2D_(std::vector< int > &vertices, std::vector< int > &geomSizes) const