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<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");
234 filenameToWrite = this->
replaceAll(filenameToWrite,
"%LEVELID",
toString(fineLevel.GetLevelID()));
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;
412 #if 0 // taw: does not compile with CGAL 4.8 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;
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;
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;
Important warning messages (one line)
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
void doAlphaHulls3D_(std::vector< int > &vertices, std::vector< int > &geomSizes) const
MueLu::DefaultLocalOrdinal LocalOrdinal
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
void writePVTU_(std::ofstream &pvtu, std::string baseFname, int numProcs) const
void writeFile_(std::ofstream &fout, std::string styleName, std::vector< int > &vertices, std::vector< int > &geomSizes) const
void doJacksPlus_(std::vector< int > &vertices, std::vector< int > &geomSizes) const
Timer to be used in factories. Similar to Monitor but with additional timers.
void doConvexHulls(std::vector< int > &vertices, std::vector< int > &geomSizes) const
One-liner description of what is happening.
Namespace for MueLu classes and methods.
void doAlphaHulls2D_(std::vector< int > &vertices, std::vector< int > &geomSizes) const
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Class that holds all level-specific information.
void doAlphaHulls_(std::vector< int > &vertices, std::vector< int > &geomSizes) const
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
Exception throws to report errors in the internal logical of the program.
void buildColormap_() const
void replaceAll(std::string &str, const std::string &from, const std::string &to)
void doGraphEdges_(std::ofstream &fout, Teuchos::RCP< Matrix > &A, Teuchos::RCP< GraphBase > &G, bool fine, int dofs) const