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" 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" 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<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");
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();
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;
195 if(doCoarseGraphEdges_)
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;
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);
210 vertex2AggId_ = vertex2AggId;
213 std::vector<GlobalOrdinal> numAggsGlobal (numProcs, 0);
214 std::vector<GlobalOrdinal> numAggsLocal (numProcs, 0);
215 std::vector<GlobalOrdinal> minGlobalAggId(numProcs, 0);
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)
221 numAggsGlobal [i] += numAggsGlobal[i-1];
222 minGlobalAggId[i] = numAggsGlobal[i-1];
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()));
238 string masterStem =
"";
241 masterStem = filenameToWrite.substr(0, filenameToWrite.rfind(
".vtu"));
242 masterStem = this->
replaceAll(masterStem,
"%PROCID",
"");
244 pvtuFilename = masterStem +
"-master.pvtu";
245 string baseFname = filenameToWrite;
247 GetOStream(
Runtime0) <<
"AggregationExportFactory: outputfile \"" << filenameToWrite <<
"\"" << std::endl;
248 ofstream fout(filenameToWrite.c_str());
249 GO numAggs = aggregates->GetNumAggregates();
252 GO indexBase = aggregates->GetMap()->getIndexBase();
253 GO offset = amalgInfo->GlobalOffset();
254 vector<GlobalOrdinal> nodeIds;
255 for (
int i = 0; i < numAggs; ++i) {
256 fout <<
"Agg " << minGlobalAggId[myRank] + i <<
" Proc " << myRank <<
":";
259 for (
int k = aggStart[i]; k < aggStart[i+1]; ++k) {
260 nodeIds.push_back((aggToRowMap[k] - offset - indexBase) / DofsPerNode + indexBase);
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());
269 for(
typename std::vector<GlobalOrdinal>::iterator printIt = nodeIds.begin(); printIt != nodeIds.end(); printIt++)
270 fout <<
" " << *printIt;
279 TEUCHOS_TEST_FOR_EXCEPTION(coords.is_null(),
Exceptions::RuntimeError,
"AggExportFactory could not get coordinates, but they are required for VTK output.");
281 numNodes_ = coords->getLocalLength();
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_)
288 cx_ = Teuchos::arcp_reinterpret_cast<
const double>(coordsCoarse->getData(0));
289 cy_ = Teuchos::arcp_reinterpret_cast<
const double>(coordsCoarse->getData(1));
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));
299 aggSizes_ = aggregates->ComputeAggregateSizes();
301 string aggStyle =
"Point Cloud";
304 aggStyle = pL.get<
string>(
"aggregation: output file: agg style");
306 catch(exception& e) {}
307 vector<int> vertices;
308 vector<int> geomSizes;
310 nodeMap_ = Amat->getMap();
311 for(LocalOrdinal i = 0; i < numNodes_; i++)
313 isRoot_.push_back(aggregates->IsRoot(i));
317 if(myRank == 0 && (numProcs != 1 || doCoarseGraphEdges_ || doFineGraphEdges_))
319 ofstream pvtu(pvtuFilename.c_str());
320 writePVTU_(pvtu, baseFname, numProcs);
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++")
328 doJacksPlus_(vertices, geomSizes);
329 else if(aggStyle ==
"Convex Hulls")
330 doConvexHulls(vertices, geomSizes);
331 else if(aggStyle ==
"Alpha Hulls")
333 #ifdef HAVE_MUELU_CGAL 334 doAlphaHulls_(vertices, geomSizes);
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);
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_);
346 writeFile_(fout, aggStyle, vertices, geomSizes);
347 if(doCoarseGraphEdges_)
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);
354 if(doFineGraphEdges_)
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);
361 if(myRank == 0 && pL.get<
bool>(
"aggregation: output file: build colormap"))
369 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
375 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
379 this->doConvexHulls2D(vertices, geomSizes, numAggs_, numNodes_, isRoot_, vertex2AggId_, xCoords_, yCoords_, zCoords_);
381 this->doConvexHulls3D(vertices, geomSizes, numAggs_, numNodes_, isRoot_, vertex2AggId_, xCoords_, yCoords_, zCoords_);
384 #ifdef HAVE_MUELU_CGAL 385 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
390 doAlphaHulls2D_(vertices, geomSizes);
392 doAlphaHulls3D_(vertices, geomSizes);
395 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
398 const double ALPHA_VAL = 2;
401 typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
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++)
414 list<Point> aggPoints;
415 vector<int> aggNodes;
416 for(
int j = 0; j < numNodes_; j++)
418 if(vertex2AggId_[j] == i)
420 Point p(xCoords_[j], yCoords_[j]);
421 aggPoints.push_back(p);
422 aggNodes.push_back(j);
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++)
432 for(
size_t k = 0; k < aggNodes.size(); k++)
434 if(fabs(segments[j][0].x == xCoords_[aggNodes[k]]) < 1e-12 && fabs(segments[j][0].y == yCoords_[aggNodes[k]]) < 1e-12)
436 vertices.push_back(aggNodes[k]);
440 for(
size_t k = 0; k < aggNodes.size(); k++)
442 if(fabs(segments[j][1].x == xCoords_[aggNodes[k]]) < 1e-12 && fabs(segments[j][1].y == yCoords_[aggNodes[k]]) < 1e-12)
444 vertices.push_back(aggNodes[k]);
448 geomSizes.push_back(2);
453 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
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;
470 for(
int i = 0; i < numAggs_; i++)
472 list<Point> aggPoints;
473 vector<int> aggNodes;
474 for(
int j = 0; j < numNodes_; j++)
476 if(vertex2AggId[j] == i)
478 Point p(xCoords_[j], yCoords_[j], zCoords_[j]);
479 aggPoints.push_back(p);
480 aggNodes.push_back(j);
483 Fixed_alpha_shape_3 hull(aggPoints.begin(), aggPoints.end(), FT(ALPHA_VAL));
484 list<Cell_handle> cells;
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++)
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++)
499 for(
size_t l = 0; l < aggNodes.size(); l++)
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)
505 vertices.push_back(aggNodes[l]);
510 geomSizes.push_back(-10);
512 for(
size_t j = 0; j < facets.size(); j++)
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]);
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();
525 for(
size_t l = 0; l < aggNodes.size(); l++)
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)
531 vertices.push_back(aggNodes[l]);
535 geomSizes.push_back(3);
537 for(
size_t j = 0; j < edges.size(); j++)
545 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
549 ArrayView<const Scalar> values;
550 ArrayView<const LocalOrdinal> neighbors;
552 vector<pair<int, int> > vert1;
553 vector<pair<int, int> > vert2;
554 if(A->isGloballyIndexed())
556 ArrayView<const GlobalOrdinal> indices;
557 for(GlobalOrdinal globRow = 0; globRow < GlobalOrdinal(A->getGlobalNumRows()); globRow++)
560 A->getGlobalRowView(globRow, indices, values);
561 neighbors = G->getNeighborVertices((LocalOrdinal) globRow);
564 while(gEdge !=
int(neighbors.size()))
568 if(neighbors[gEdge] == indices[aEdge])
571 vert1.push_back(pair<int, int>(
int(globRow), neighbors[gEdge]));
578 vert2.push_back(pair<int, int>(
int(globRow), neighbors[gEdge]));
584 vert1.push_back(pair<int, int>(
int(globRow), neighbors[gEdge]));
592 ArrayView<const LocalOrdinal> indices;
593 for(LocalOrdinal locRow = 0; locRow < LocalOrdinal(A->getNodeNumRows()); locRow++)
596 A->getLocalRowView(locRow, indices, values);
597 neighbors = G->getNeighborVertices(locRow);
601 while(gEdge !=
int(neighbors.size()))
605 if(neighbors[gEdge] == indices[aEdge])
607 vert1.push_back(pair<int, int>(locRow, neighbors[gEdge]));
613 vert2.push_back(pair<int, int>(locRow, neighbors[gEdge]));
619 vert1.push_back(pair<int, int>(locRow, neighbors[gEdge]));
625 for(
size_t i = 0; i < vert1.size(); i ++)
627 if(vert1[i].first > vert1[i].second)
629 int temp = vert1[i].first;
630 vert1[i].first = vert1[i].second;
631 vert1[i].second = temp;
634 for(
size_t i = 0; i < vert2.size(); i++)
636 if(vert2[i].first > vert2[i].second)
638 int temp = vert2[i].first;
639 vert2[i].first = vert2[i].second;
640 vert2[i].second = temp;
643 sort(vert1.begin(), vert1.end());
644 vector<pair<int, int> >::iterator newEnd = unique(vert1.begin(), vert1.end());
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());
650 points1.reserve(2 * vert1.size());
651 for(
size_t i = 0; i < vert1.size(); i++)
653 points1.push_back(vert1[i].first);
654 points1.push_back(vert1[i].second);
657 points2.reserve(2 * vert2.size());
658 for(
size_t i = 0; i < vert2.size(); i++)
660 points2.push_back(vert2[i].first);
661 points2.push_back(vert2[i].second);
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;
672 for(
size_t i = 0; i < unique1.size(); i++)
674 fout << CONTRAST_1_ <<
" ";
676 fout << endl << indent;
678 for(
size_t i = 0; i < unique2.size(); i++)
680 fout << CONTRAST_2_ <<
" ";
681 if((i + 2 * vert1.size()) % 25 == 24)
682 fout << endl << indent;
685 fout <<
" </DataArray>" << endl;
686 fout <<
" <DataArray type=\"Int32\" Name=\"Aggregate\" format=\"ascii\">" << endl;
688 for(
size_t i = 0; i < unique1.size(); i++)
690 fout << CONTRAST_1_ <<
" ";
692 fout << endl << indent;
694 for(
size_t i = 0; i < unique2.size(); i++)
696 fout << CONTRAST_2_ <<
" ";
697 if((i + 2 * vert2.size()) % 25 == 24)
698 fout << endl << indent;
701 fout <<
" </DataArray>" << endl;
702 fout <<
" <DataArray type=\"Int32\" Name=\"Processor\" format=\"ascii\">" << endl;
704 for(
size_t i = 0; i < unique1.size() + unique2.size(); i++)
706 fout << myRank_ <<
" ";
708 fout << endl << indent;
711 fout <<
" </DataArray>" << endl;
712 fout <<
" </PointData>" << endl;
713 fout <<
" <Points>" << endl;
714 fout <<
" <DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"ascii\">" << endl;
716 for(
size_t i = 0; i < unique1.size(); i++)
720 fout << xCoords_[unique1[i]] <<
" " << yCoords_[unique1[i]] <<
" ";
722 fout << zCoords_[unique1[i]] <<
" ";
726 fout << endl << indent;
730 fout << cx_[unique1[i]] <<
" " << cy_[unique1[i]] <<
" ";
732 fout << cz_[unique1[i]] <<
" ";
736 fout << endl << indent;
739 for(
size_t i = 0; i < unique2.size(); i++)
743 fout << xCoords_[unique2[i]] <<
" " << yCoords_[unique2[i]] <<
" ";
745 fout << zCoords_[unique2[i]] <<
" ";
749 fout << endl << indent;
753 fout << cx_[unique2[i]] <<
" " << cy_[unique2[i]] <<
" ";
755 fout << cz_[unique2[i]] <<
" ";
758 if((i + unique1.size()) % 2)
759 fout << endl << indent;
763 fout <<
" </DataArray>" << endl;
764 fout <<
" </Points>" << endl;
765 fout <<
" <Cells>" << endl;
766 fout <<
" <DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">" << endl;
768 for(
size_t i = 0; i < points1.size(); i++)
770 fout << points1[i] <<
" ";
772 fout << endl << indent;
774 for(
size_t i = 0; i < points2.size(); i++)
776 fout << points2[i] + unique1.size() <<
" ";
777 if((i + points1.size()) % 10 == 9)
778 fout << endl << indent;
781 fout <<
" </DataArray>" << endl;
782 fout <<
" <DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">" << endl;
785 for(
size_t i = 0; i < vert1.size() + vert2.size(); i++)
788 fout << offset <<
" ";
790 fout << endl << indent;
793 fout <<
" </DataArray>" << endl;
794 fout <<
" <DataArray type=\"Int32\" Name=\"types\" format=\"ascii\">" << endl;
796 for(
size_t i = 0; i < vert1.size() + vert2.size(); i++)
800 fout << endl << indent;
803 fout <<
" </DataArray>" << endl;
804 fout <<
" </Cells>" << endl;
805 fout <<
" </Piece>" << endl;
806 fout <<
" </UnstructuredGrid>" << endl;
807 fout <<
"</VTKFile>" << endl;
810 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
814 vector<int> uniqueFine = this->makeUnique(vertices);
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;
824 bool localIsGlobal = GlobalOrdinal(nodeMap_->getGlobalNumElements()) == GlobalOrdinal(nodeMap_->getNodeNumElements());
825 for(
size_t i = 0; i < uniqueFine.size(); i++)
829 fout << uniqueFine[i] <<
" ";
832 fout << nodeMap_->getGlobalElement(uniqueFine[i]) <<
" ";
834 fout << endl << indent;
837 fout <<
" </DataArray>" << endl;
838 fout <<
" <DataArray type=\"Int32\" Name=\"Aggregate\" format=\"ascii\">" << endl;
840 for(
size_t i = 0; i < uniqueFine.size(); i++)
842 fout << aggsOffset_ + vertex2AggId_[uniqueFine[i]] <<
" ";
844 fout << endl << indent;
847 fout <<
" </DataArray>" << endl;
848 fout <<
" <DataArray type=\"Int32\" Name=\"Processor\" format=\"ascii\">" << endl;
850 for(
size_t i = 0; i < uniqueFine.size(); i++)
852 fout << myRank_ <<
" ";
854 fout << endl << indent;
857 fout <<
" </DataArray>" << endl;
858 fout <<
" </PointData>" << endl;
859 fout <<
" <Points>" << endl;
860 fout <<
" <DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"ascii\">" << endl;
862 for(
size_t i = 0; i < uniqueFine.size(); i++)
864 fout << xCoords_[uniqueFine[i]] <<
" " << yCoords_[uniqueFine[i]] <<
" ";
868 fout << zCoords_[uniqueFine[i]] <<
" ";
870 fout << endl << indent;
873 fout <<
" </DataArray>" << endl;
874 fout <<
" </Points>" << endl;
875 fout <<
" <Cells>" << endl;
876 fout <<
" <DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">" << endl;
878 for(
size_t i = 0; i < vertices.size(); i++)
880 fout << vertices[i] <<
" ";
882 fout << endl << indent;
885 fout <<
" </DataArray>" << endl;
886 fout <<
" <DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">" << endl;
889 for(
size_t i = 0; i < geomSizes.size(); i++)
891 accum += geomSizes[i];
892 fout << accum <<
" ";
894 fout << endl << indent;
897 fout <<
" </DataArray>" << endl;
898 fout <<
" <DataArray type=\"Int32\" Name=\"types\" format=\"ascii\">" << endl;
900 for(
size_t i = 0; i < geomSizes.size(); i++)
917 fout << endl << indent;
920 fout <<
" </DataArray>" << endl;
921 fout <<
" </Cells>" << endl;
922 fout <<
" </Piece>" << endl;
923 fout <<
" </UnstructuredGrid>" << endl;
924 fout <<
"</VTKFile>" << endl;
927 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
933 ofstream color(
"random-colormap.xml");
934 color <<
"<ColorMap name=\"MueLu-Random\" space=\"RGB\">" << endl;
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;
941 for(
int i = 0; i < 5000; i += 4)
943 color <<
" <Point x=\"" << i <<
"\" o=\"1\" r=\"" << (rand() % 50) / 256.0 <<
"\" g=\"" << (rand() % 256) / 256.0 <<
"\" b=\"" << (rand() % 256) / 256.0 <<
"\"/>" << endl;
945 color <<
"</ColorMap>" << endl;
950 GetOStream(
Warnings0) <<
" Error while building colormap file: " << e.what() << endl;
954 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
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++)
974 pvtu <<
" <Piece Source=\"" << this->
replaceAll(baseFname,
"%PROCID",
toString(i)) <<
"\"/>" << endl;
977 if(doFineGraphEdges_)
979 for(
int i = 0; i < numProcs; i++)
982 pvtu <<
" <Piece Source=\"" << fn.insert(fn.rfind(
".vtu"),
"-finegraph") <<
"\"/>" << endl;
985 if(doCoarseGraphEdges_)
987 for(
int i = 0; i < numProcs; i++)
990 pvtu <<
" <Piece Source=\"" << fn.insert(fn.rfind(
".vtu"),
"-coarsegraph") <<
"\"/>" << endl;
993 pvtu <<
" </PUnstructuredGrid>" << endl;
994 pvtu <<
"</VTKFile>" << endl;
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.
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
void buildColormap_() const