MueLu Version of the Day
Loading...
Searching...
No Matches
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
85namespace 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->getLocalNumRows()); 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_->getLocalNumElements());
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_ */
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultGlobalOrdinal GlobalOrdinal
MueLu::DefaultNode Node
void doJacksPlus_(std::vector< int > &vertices, std::vector< int > &geomSizes) const
void writePVTU_(std::ofstream &pvtu, std::string baseFname, int numProcs) const
void doConvexHulls(std::vector< int > &vertices, std::vector< int > &geomSizes) const
void doGraphEdges_(std::ofstream &fout, Teuchos::RCP< Matrix > &A, Teuchos::RCP< GraphBase > &G, bool fine, int dofs) const
void doAlphaHulls_(std::vector< int > &vertices, std::vector< int > &geomSizes) const
void writeFile_(std::ofstream &fout, std::string styleName, std::vector< int > &vertices, std::vector< int > &geomSizes) const
void doAlphaHulls2D_(std::vector< int > &vertices, std::vector< int > &geomSizes) const
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
void doAlphaHulls3D_(std::vector< int > &vertices, std::vector< int > &geomSizes) const
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:76
Namespace for MueLu classes and methods.
@ Warnings0
Important warning messages (one line)
@ Runtime0
One-liner description of what is happening.
void replaceAll(std::string &str, const std::string &from, const std::string &to)
std::string toString(const T &what)
Little helper function to convert non-string types to strings.