53#ifndef MUELU_AGGREGATIONEXPORTFACTORY_DEF_HPP_
54#define MUELU_AGGREGATIONEXPORTFACTORY_DEF_HPP_
56#include <Xpetra_Matrix.hpp>
57#include <Xpetra_CrsMatrixWrap.hpp>
58#include <Xpetra_ImportFactory.hpp>
59#include <Xpetra_MultiVectorFactory.hpp>
62#include "MueLu_Aggregates.hpp"
63#include "MueLu_Graph.hpp"
64#include "MueLu_AmalgamationFactory.hpp"
65#include "MueLu_AmalgamationInfo.hpp"
67#include "MueLu_Utilities.hpp"
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"
87 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
89 RCP<ParameterList> validParamList = rcp(
new ParameterList());
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";
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.");
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");
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");
109 validParamList->set<
int > (
"aggregation: output file: iter", 0,
"nonlinear iteration variable for output file name");
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;
117 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
119 Input(fineLevel,
"Aggregates");
120 Input(fineLevel,
"DofsPerNode");
121 Input(fineLevel,
"UnAmalgamationInfo");
123 const ParameterList & pL = GetParameterList();
125 if(pL.isParameter(
"aggregation: output filename") && pL.get<std::string>(
"aggregation: output filename").length())
127 Input(fineLevel,
"Coordinates");
128 Input(fineLevel,
"A");
129 Input(fineLevel,
"Graph");
130 if(pL.get<
bool>(
"aggregation: output file: coarse graph edges"))
132 Input(coarseLevel,
"Coordinates");
133 Input(coarseLevel,
"A");
134 Input(coarseLevel,
"Graph");
139 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
143 const ParameterList& pL = GetParameterList();
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");
150 string pvtuFilename =
"";
151 string localFilename = pL.get<std::string>(
"Output filename");
152 string filenameToWrite;
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())
159 filenameToWrite = masterFilename;
160 if(filenameToWrite.rfind(
".vtu") == string::npos)
161 filenameToWrite.append(
".vtu");
162 if(numProcs > 1 && filenameToWrite.rfind(
"%PROCID") == string::npos)
163 filenameToWrite.insert(filenameToWrite.rfind(
".vtu"),
"-proc%PROCID");
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");
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();
189 if (aggregates->AggregatesCrossProcessors())
191 RCP<Import> coordImporter = Xpetra::ImportFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(coords->getMap(), Amat->getColMap());
193 ghostedCoords->doImport(*coords, *coordImporter, Xpetra::INSERT);
194 coords = ghostedCoords;
196 if(doCoarseGraphEdges_)
198 RCP<Import> coordImporter = Xpetra::ImportFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(coordsCoarse->getMap(), Ac->getColMap());
200 ghostedCoords->doImport(*coordsCoarse, *coordImporter, Xpetra::INSERT);
201 coordsCoarse = ghostedCoords;
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);
211 vertex2AggId_ = vertex2AggId;
214 std::vector<GlobalOrdinal> numAggsGlobal (numProcs, 0);
215 std::vector<GlobalOrdinal> numAggsLocal (numProcs, 0);
216 std::vector<GlobalOrdinal> minGlobalAggId(numProcs, 0);
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)
222 numAggsGlobal [i] += numAggsGlobal[i-1];
223 minGlobalAggId[i] = numAggsGlobal[i-1];
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");
239 string masterStem =
"";
242 masterStem = filenameToWrite.substr(0, filenameToWrite.rfind(
".vtu"));
243 masterStem = this->
replaceAll(masterStem,
"%PROCID",
"");
245 pvtuFilename = masterStem +
"-master.pvtu";
246 string baseFname = filenameToWrite;
248 GetOStream(
Runtime0) <<
"AggregationExportFactory: outputfile \"" << filenameToWrite <<
"\"" << std::endl;
249 ofstream fout(filenameToWrite.c_str());
250 GO numAggs = aggregates->GetNumAggregates();
253 GO indexBase = aggregates->GetMap()->getIndexBase();
254 GO offset = amalgInfo->GlobalOffset();
255 vector<GlobalOrdinal> nodeIds;
256 for (
int i = 0; i < numAggs; ++i) {
257 fout <<
"Agg " << minGlobalAggId[myRank] + i <<
" Proc " << myRank <<
":";
260 for (
int k = aggStart[i]; k < aggStart[i+1]; ++k) {
261 nodeIds.push_back((aggToRowMap[k] - offset - indexBase) / DofsPerNode + indexBase);
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());
270 for(
typename std::vector<GlobalOrdinal>::iterator printIt = nodeIds.begin(); printIt != nodeIds.end(); printIt++)
271 fout <<
" " << *printIt;
280 TEUCHOS_TEST_FOR_EXCEPTION(coords.is_null(),
Exceptions::RuntimeError,
"AggExportFactory could not get coordinates, but they are required for VTK output.");
282 numNodes_ = coords->getLocalLength();
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_)
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));
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));
300 aggSizes_ = aggregates->ComputeAggregateSizes();
302 string aggStyle =
"Point Cloud";
305 aggStyle = pL.get<
string>(
"aggregation: output file: agg style");
307 catch(std::exception& e) {}
308 vector<int> vertices;
309 vector<int> geomSizes;
311 nodeMap_ = Amat->getMap();
314 isRoot_.push_back(aggregates->IsRoot(i));
318 if(myRank == 0 && (numProcs != 1 || doCoarseGraphEdges_ || doFineGraphEdges_))
320 ofstream pvtu(pvtuFilename.c_str());
321 writePVTU_(pvtu, baseFname, numProcs);
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++")
329 doJacksPlus_(vertices, geomSizes);
330 else if(aggStyle ==
"Convex Hulls")
331 doConvexHulls(vertices, geomSizes);
332 else if(aggStyle ==
"Alpha Hulls")
334 #ifdef HAVE_MUELU_CGAL
335 doAlphaHulls_(vertices, geomSizes);
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);
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_);
347 writeFile_(fout, aggStyle, vertices, geomSizes);
348 if(doCoarseGraphEdges_)
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);
355 if(doFineGraphEdges_)
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);
362 if(myRank == 0 && pL.get<
bool>(
"aggregation: output file: build colormap"))
370 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
376 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
380 this->doConvexHulls2D(vertices, geomSizes, numAggs_, numNodes_, isRoot_, vertex2AggId_, xCoords_, yCoords_);
382 this->doConvexHulls3D(vertices, geomSizes, numAggs_, numNodes_, isRoot_, vertex2AggId_, xCoords_, yCoords_, zCoords_);
385#ifdef HAVE_MUELU_CGAL
386 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
391 doAlphaHulls2D_(vertices, geomSizes);
393 doAlphaHulls3D_(vertices, geomSizes);
396 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
402 typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
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;
413 for(
int i = 0; i < numAggs_; i++)
416 list<Point> aggPoints;
417 vector<int> aggNodes;
418 for(
int j = 0; j < numNodes_; j++)
420 if(vertex2AggId_[j] == i)
422 Point p(xCoords_[j], yCoords_[j]);
423 aggPoints.push_back(p);
424 aggNodes.push_back(j);
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++)
434 for(
size_t k = 0; k < aggNodes.size(); k++)
436 if(fabs(segments[j][0].x == xCoords_[aggNodes[k]]) < 1e-12 && fabs(segments[j][0].y == yCoords_[aggNodes[k]]) < 1e-12)
438 vertices.push_back(aggNodes[k]);
442 for(
size_t k = 0; k < aggNodes.size(); k++)
444 if(fabs(segments[j][1].x == xCoords_[aggNodes[k]]) < 1e-12 && fabs(segments[j][1].y == yCoords_[aggNodes[k]]) < 1e-12)
446 vertices.push_back(aggNodes[k]);
450 geomSizes.push_back(2);
456 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
459 typedef CGAL::Exact_predicates_inexact_constructions_kernel Gt;
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;
475 for(
int i = 0; i < numAggs_; i++)
477 list<Point> aggPoints;
478 vector<int> aggNodes;
479 for(
int j = 0; j < numNodes_; j++)
481 if(vertex2AggId[j] == i)
483 Point p(xCoords_[j], yCoords_[j], zCoords_[j]);
484 aggPoints.push_back(p);
485 aggNodes.push_back(j);
488 Fixed_alpha_shape_3 hull(aggPoints.begin(), aggPoints.end(), FT(ALPHA_VAL));
489 list<Cell_handle> cells;
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++)
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++)
504 for(
size_t l = 0; l < aggNodes.size(); l++)
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)
510 vertices.push_back(aggNodes[l]);
515 geomSizes.push_back(-10);
517 for(
size_t j = 0; j < facets.size(); j++)
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]);
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();
530 for(
size_t l = 0; l < aggNodes.size(); l++)
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)
536 vertices.push_back(aggNodes[l]);
540 geomSizes.push_back(3);
542 for(
size_t j = 0; j < edges.size(); j++)
551 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
555 ArrayView<const Scalar> values;
556 ArrayView<const LocalOrdinal> neighbors;
558 vector<pair<int, int> > vert1;
559 vector<pair<int, int> > vert2;
560 if(A->isGloballyIndexed())
562 ArrayView<const GlobalOrdinal> indices;
566 A->getGlobalRowView(globRow, indices, values);
567 neighbors = G->getNeighborVertices((
LocalOrdinal) globRow);
570 while(gEdge !=
int(neighbors.size()))
574 if(neighbors[gEdge] == indices[aEdge])
577 vert1.push_back(pair<int, int>(
int(globRow), neighbors[gEdge]));
584 vert2.push_back(pair<int, int>(
int(globRow), neighbors[gEdge]));
590 vert1.push_back(pair<int, int>(
int(globRow), neighbors[gEdge]));
598 ArrayView<const LocalOrdinal> indices;
602 A->getLocalRowView(locRow, indices, values);
603 neighbors = G->getNeighborVertices(locRow);
607 while(gEdge !=
int(neighbors.size()))
611 if(neighbors[gEdge] == indices[aEdge])
613 vert1.push_back(pair<int, int>(locRow, neighbors[gEdge]));
619 vert2.push_back(pair<int, int>(locRow, neighbors[gEdge]));
625 vert1.push_back(pair<int, int>(locRow, neighbors[gEdge]));
631 for(
size_t i = 0; i < vert1.size(); i ++)
633 if(vert1[i].first > vert1[i].second)
635 int temp = vert1[i].first;
636 vert1[i].first = vert1[i].second;
637 vert1[i].second = temp;
640 for(
size_t i = 0; i < vert2.size(); i++)
642 if(vert2[i].first > vert2[i].second)
644 int temp = vert2[i].first;
645 vert2[i].first = vert2[i].second;
646 vert2[i].second = temp;
649 sort(vert1.begin(), vert1.end());
650 vector<pair<int, int> >::iterator newEnd = unique(vert1.begin(), vert1.end());
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());
656 points1.reserve(2 * vert1.size());
657 for(
size_t i = 0; i < vert1.size(); i++)
659 points1.push_back(vert1[i].first);
660 points1.push_back(vert1[i].second);
663 points2.reserve(2 * vert2.size());
664 for(
size_t i = 0; i < vert2.size(); i++)
666 points2.push_back(vert2[i].first);
667 points2.push_back(vert2[i].second);
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;
678 for(
size_t i = 0; i < unique1.size(); i++)
680 fout << CONTRAST_1_ <<
" ";
682 fout << endl << indent;
684 for(
size_t i = 0; i < unique2.size(); i++)
686 fout << CONTRAST_2_ <<
" ";
687 if((i + 2 * vert1.size()) % 25 == 24)
688 fout << endl << indent;
691 fout <<
" </DataArray>" << endl;
692 fout <<
" <DataArray type=\"Int32\" Name=\"Aggregate\" format=\"ascii\">" << endl;
694 for(
size_t i = 0; i < unique1.size(); i++)
696 fout << CONTRAST_1_ <<
" ";
698 fout << endl << indent;
700 for(
size_t i = 0; i < unique2.size(); i++)
702 fout << CONTRAST_2_ <<
" ";
703 if((i + 2 * vert2.size()) % 25 == 24)
704 fout << endl << indent;
707 fout <<
" </DataArray>" << endl;
708 fout <<
" <DataArray type=\"Int32\" Name=\"Processor\" format=\"ascii\">" << endl;
710 for(
size_t i = 0; i < unique1.size() + unique2.size(); i++)
712 fout << myRank_ <<
" ";
714 fout << endl << indent;
717 fout <<
" </DataArray>" << endl;
718 fout <<
" </PointData>" << endl;
719 fout <<
" <Points>" << endl;
720 fout <<
" <DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"ascii\">" << endl;
722 for(
size_t i = 0; i < unique1.size(); i++)
726 fout << xCoords_[unique1[i]] <<
" " << yCoords_[unique1[i]] <<
" ";
728 fout << zCoords_[unique1[i]] <<
" ";
732 fout << endl << indent;
736 fout << cx_[unique1[i]] <<
" " << cy_[unique1[i]] <<
" ";
738 fout << cz_[unique1[i]] <<
" ";
742 fout << endl << indent;
745 for(
size_t i = 0; i < unique2.size(); i++)
749 fout << xCoords_[unique2[i]] <<
" " << yCoords_[unique2[i]] <<
" ";
751 fout << zCoords_[unique2[i]] <<
" ";
755 fout << endl << indent;
759 fout << cx_[unique2[i]] <<
" " << cy_[unique2[i]] <<
" ";
761 fout << cz_[unique2[i]] <<
" ";
764 if((i + unique1.size()) % 2)
765 fout << endl << indent;
769 fout <<
" </DataArray>" << endl;
770 fout <<
" </Points>" << endl;
771 fout <<
" <Cells>" << endl;
772 fout <<
" <DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">" << endl;
774 for(
size_t i = 0; i < points1.size(); i++)
776 fout << points1[i] <<
" ";
778 fout << endl << indent;
780 for(
size_t i = 0; i < points2.size(); i++)
782 fout << points2[i] + unique1.size() <<
" ";
783 if((i + points1.size()) % 10 == 9)
784 fout << endl << indent;
787 fout <<
" </DataArray>" << endl;
788 fout <<
" <DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">" << endl;
791 for(
size_t i = 0; i < vert1.size() + vert2.size(); i++)
794 fout << offset <<
" ";
796 fout << endl << indent;
799 fout <<
" </DataArray>" << endl;
800 fout <<
" <DataArray type=\"Int32\" Name=\"types\" format=\"ascii\">" << endl;
802 for(
size_t i = 0; i < vert1.size() + vert2.size(); i++)
806 fout << endl << indent;
809 fout <<
" </DataArray>" << endl;
810 fout <<
" </Cells>" << endl;
811 fout <<
" </Piece>" << endl;
812 fout <<
" </UnstructuredGrid>" << endl;
813 fout <<
"</VTKFile>" << endl;
816 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
820 vector<int> uniqueFine = this->makeUnique(vertices);
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;
831 for(
size_t i = 0; i < uniqueFine.size(); i++)
835 fout << uniqueFine[i] <<
" ";
838 fout << nodeMap_->getGlobalElement(uniqueFine[i]) <<
" ";
840 fout << endl << indent;
843 fout <<
" </DataArray>" << endl;
844 fout <<
" <DataArray type=\"Int32\" Name=\"Aggregate\" format=\"ascii\">" << endl;
846 for(
size_t i = 0; i < uniqueFine.size(); i++)
848 if(vertex2AggId_[uniqueFine[i]]==-1)
849 fout << vertex2AggId_[uniqueFine[i]] <<
" ";
851 fout << aggsOffset_ + vertex2AggId_[uniqueFine[i]] <<
" ";
853 fout << endl << indent;
856 fout <<
" </DataArray>" << endl;
857 fout <<
" <DataArray type=\"Int32\" Name=\"Processor\" format=\"ascii\">" << endl;
859 for(
size_t i = 0; i < uniqueFine.size(); i++)
861 fout << myRank_ <<
" ";
863 fout << endl << indent;
866 fout <<
" </DataArray>" << endl;
867 fout <<
" </PointData>" << endl;
868 fout <<
" <Points>" << endl;
869 fout <<
" <DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"ascii\">" << endl;
871 for(
size_t i = 0; i < uniqueFine.size(); i++)
873 fout << xCoords_[uniqueFine[i]] <<
" " << yCoords_[uniqueFine[i]] <<
" ";
877 fout << zCoords_[uniqueFine[i]] <<
" ";
879 fout << endl << indent;
882 fout <<
" </DataArray>" << endl;
883 fout <<
" </Points>" << endl;
884 fout <<
" <Cells>" << endl;
885 fout <<
" <DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">" << endl;
887 for(
size_t i = 0; i < vertices.size(); i++)
889 fout << vertices[i] <<
" ";
891 fout << endl << indent;
894 fout <<
" </DataArray>" << endl;
895 fout <<
" <DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">" << endl;
898 for(
size_t i = 0; i < geomSizes.size(); i++)
900 accum += geomSizes[i];
901 fout << accum <<
" ";
903 fout << endl << indent;
906 fout <<
" </DataArray>" << endl;
907 fout <<
" <DataArray type=\"Int32\" Name=\"types\" format=\"ascii\">" << endl;
909 for(
size_t i = 0; i < geomSizes.size(); i++)
926 fout << endl << indent;
929 fout <<
" </DataArray>" << endl;
930 fout <<
" </Cells>" << endl;
931 fout <<
" </Piece>" << endl;
932 fout <<
" </UnstructuredGrid>" << endl;
933 fout <<
"</VTKFile>" << endl;
936 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
942 ofstream color(
"random-colormap.xml");
943 color <<
"<ColorMap name=\"MueLu-Random\" space=\"RGB\">" << endl;
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;
950 for(
int i = 0; i < 5000; i += 4)
952 color <<
" <Point x=\"" << i <<
"\" o=\"1\" r=\"" << (rand() % 50) / 256.0 <<
"\" g=\"" << (rand() % 256) / 256.0 <<
"\" b=\"" << (rand() % 256) / 256.0 <<
"\"/>" << endl;
954 color <<
"</ColorMap>" << endl;
957 catch(std::exception& e)
959 GetOStream(
Warnings0) <<
" Error while building colormap file: " << e.what() << endl;
963 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
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++)
983 pvtu <<
" <Piece Source=\"" << this->
replaceAll(baseFname,
"%PROCID",
toString(i)) <<
"\"/>" << endl;
986 if(doFineGraphEdges_)
988 for(
int i = 0; i < numProcs; i++)
991 pvtu <<
" <Piece Source=\"" << fn.insert(fn.rfind(
".vtu"),
"-finegraph") <<
"\"/>" << endl;
994 if(doCoarseGraphEdges_)
996 for(
int i = 0; i < numProcs; i++)
999 pvtu <<
" <Piece Source=\"" << fn.insert(fn.rfind(
".vtu"),
"-coarsegraph") <<
"\"/>" << endl;
1002 pvtu <<
" </PUnstructuredGrid>" << endl;
1003 pvtu <<
"</VTKFile>" << endl;
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultGlobalOrdinal GlobalOrdinal
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 buildColormap_() 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.
int GetLevelID() const
Return level number.
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.