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"))
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);
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++)
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;