MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_AggregationExportFactory_def.hpp
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// MueLu: A package for multigrid based preconditioning
4//
5// Copyright 2012 NTESS and the MueLu contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
10/*
11 * MueLu_AggregationExportFactory_def.hpp
12 *
13 * Created on: Feb 10, 2012
14 * Author: wiesner
15 */
16
17#ifndef MUELU_AGGREGATIONEXPORTFACTORY_DEF_HPP_
18#define MUELU_AGGREGATIONEXPORTFACTORY_DEF_HPP_
19
20#include <Xpetra_Matrix.hpp>
21#include <Xpetra_CrsMatrixWrap.hpp>
22#include <Xpetra_ImportFactory.hpp>
23#include <Xpetra_MultiVectorFactory.hpp>
24#include "MueLu_AggregationExportFactory_decl.hpp"
25#include "MueLu_Level.hpp"
26#include "MueLu_Aggregates.hpp"
27#include "MueLu_FactoryManagerBase.hpp"
28
29#include "MueLu_AmalgamationFactory.hpp"
30#include "MueLu_AmalgamationInfo.hpp"
31#include "MueLu_Monitor.hpp"
32#include <vector>
33#include <list>
34#include <algorithm>
35#include <string>
36#include <stdexcept>
37#include <cstdio>
38#include <cmath>
39
40namespace MueLu {
41
42template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
43AggregationExportFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::AggregationExportFactory()
44 : doFineGraphEdges_(false)
45 , doCoarseGraphEdges_(false)
46 , numNodes_(0)
47 , numAggs_(0)
48 , dims_(0)
49 , myRank_(-1)
50 , aggsOffset_(0) {}
51
52template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
53AggregationExportFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::~AggregationExportFactory() = default;
54
55template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
56RCP<const ParameterList> AggregationExportFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetValidParameterList() const {
57 RCP<ParameterList> validParamList = rcp(new ParameterList());
58
59 std::string output_msg =
60 "Output filename template (%TIMESTEP is replaced by \'Output file: time step\' variable,"
61 "%ITER is replaced by \'Output file: iter\' variable, %LEVELID is replaced level id, %PROCID is replaced by processor id)";
62 std::string output_def = "aggs_level%LEVELID_proc%PROCID.out";
63
64 validParamList->set<RCP<const FactoryBase>>("A", Teuchos::null, "Factory for A.");
65 validParamList->set<RCP<const FactoryBase>>("Coordinates", Teuchos::null, "Factory for Coordinates.");
66 validParamList->set<RCP<const FactoryBase>>("Material", Teuchos::null, "Factory for Material.");
67 validParamList->set<RCP<const FactoryBase>>("Graph", Teuchos::null, "Factory for Graph.");
68 validParamList->set<RCP<const FactoryBase>>("Aggregates", Teuchos::null, "Factory for Aggregates.");
69 validParamList->set<RCP<const FactoryBase>>("AggregateQualities", Teuchos::null, "Factory for AggregateQualities.");
70 validParamList->set<RCP<const FactoryBase>>("UnAmalgamationInfo", Teuchos::null, "Factory for UnAmalgamationInfo.");
71 validParamList->set<RCP<const FactoryBase>>("DofsPerNode", Teuchos::null, "Factory for DofsPerNode.");
72 // CMS/BMK: Old style factory-only options. Deprecate me.
73 validParamList->set<std::string>("Output filename", output_def, output_msg);
74 validParamList->set<int>("Output file: time step", 0, "time step variable for output file name");
75 validParamList->set<int>("Output file: iter", 0, "nonlinear iteration variable for output file name");
76
77 // New-style master list options (here are same defaults as in masterList.xml)
78 validParamList->set<std::string>("aggregation: output filename", "", "filename for VTK-style visualization output");
79 validParamList->set<int>("aggregation: output file: time step", 0, "time step variable for output file name"); // Remove me?
80 validParamList->set<int>("aggregation: output file: iter", 0, "nonlinear iteration variable for output file name"); // Remove me?
81 validParamList->set<std::string>("aggregation: output file: agg style", "Point Cloud", "style of aggregate visualization for VTK output");
82 validParamList->set<bool>("aggregation: output file: fine graph edges", false, "Whether to draw all fine node connections along with the aggregates.");
83 validParamList->set<bool>("aggregation: output file: coarse graph edges", false, "Whether to draw all coarse node connections along with the aggregates.");
84 validParamList->set<bool>("aggregation: output file: build colormap", false, "Whether to output a random colormap for ParaView in a separate XML file.");
85 validParamList->set<bool>("aggregation: output file: aggregate qualities", false, "Whether to plot the aggregate quality.");
86 validParamList->set<bool>("aggregation: output file: material", false, "Whether to plot the material.");
87 return validParamList;
88}
89
90template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
91void AggregationExportFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::DeclareInput(Level& fineLevel, Level& coarseLevel) const {
92 Input(fineLevel, "Aggregates"); //< factory which created aggregates
93 Input(fineLevel, "DofsPerNode"); //< CoalesceAndDropFactory (needed for DofsPerNode variable)
94 Input(fineLevel, "UnAmalgamationInfo"); //< AmalgamationFactory (needed for UnAmalgamationInfo variable)
95
96 const ParameterList& pL = GetParameterList();
97 // Only pull in coordinates if the user explicitly requests direct VTK output, so as not to break uses of old code
98 if (pL.isParameter("aggregation: output filename") && pL.get<std::string>("aggregation: output filename").length()) {
99 Input(fineLevel, "Coordinates");
100 Input(fineLevel, "A");
101 Input(fineLevel, "Graph");
102 if (pL.get<bool>("aggregation: output file: coarse graph edges")) {
103 Input(coarseLevel, "Coordinates");
104 Input(coarseLevel, "A");
105 Input(coarseLevel, "Graph");
106 }
107 }
108
109 if (pL.get<bool>("aggregation: output file: aggregate qualities")) {
110 Input(coarseLevel, "AggregateQualities");
111 }
112
113 if (pL.get<bool>("aggregation: output file: material")) {
114 Input(fineLevel, "Material");
115 }
116}
117
118template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
119void AggregationExportFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(Level& fineLevel, Level& coarseLevel) const {
120 using namespace std;
121 // Decide which build function to follow, based on input params
122 const ParameterList& pL = GetParameterList();
123 FactoryMonitor m(*this, "AggregationExportFactory", coarseLevel);
124 Teuchos::RCP<Aggregates> aggregates = Get<Teuchos::RCP<Aggregates>>(fineLevel, "Aggregates");
125 Teuchos::RCP<const Teuchos::Comm<int>> comm = aggregates->GetMap()->getComm();
126 int numProcs = comm->getSize();
127 int myRank = comm->getRank();
128 string masterFilename = pL.get<std::string>("aggregation: output filename"); // filename parameter from master list
129 string pvtuFilename = ""; // only root processor will set this
130 string localFilename = pL.get<std::string>("Output filename");
131 string filenameToWrite;
132 bool useVTK = false;
133 doCoarseGraphEdges_ = pL.get<bool>("aggregation: output file: coarse graph edges");
134 doFineGraphEdges_ = pL.get<bool>("aggregation: output file: fine graph edges");
135 doAggQuality_ = pL.get<bool>("aggregation: output file: aggregate qualities");
136 doMaterial_ = pL.get<bool>("aggregation: output file: material");
137 if (masterFilename.length()) {
138 useVTK = true;
139 filenameToWrite = masterFilename;
140 if (filenameToWrite.rfind(".vtu") == string::npos) // Must have the file extension in the name
141 filenameToWrite.append(".vtu");
142 if (numProcs > 1 && filenameToWrite.rfind("%PROCID") == string::npos) // filename can't be identical between processsors in parallel problem
143 filenameToWrite.insert(filenameToWrite.rfind(".vtu"), "-proc%PROCID");
144 } else
145 filenameToWrite = localFilename;
146 LocalOrdinal DofsPerNode = Get<LocalOrdinal>(fineLevel, "DofsPerNode");
147 Teuchos::RCP<AmalgamationInfo> amalgInfo = Get<RCP<AmalgamationInfo>>(fineLevel, "UnAmalgamationInfo");
148 Teuchos::RCP<Matrix> Amat = Get<RCP<Matrix>>(fineLevel, "A");
149 Teuchos::RCP<Matrix> Ac;
150 if (doCoarseGraphEdges_)
151 Ac = Get<RCP<Matrix>>(coarseLevel, "A");
152 Teuchos::RCP<CoordinateMultiVector> coords = Teuchos::null;
153 Teuchos::RCP<CoordinateMultiVector> coordsCoarse = Teuchos::null;
154 if (doAggQuality_)
155 qualities_ = Get<Teuchos::RCP<MultiVector>>(coarseLevel, "AggregateQualities");
156 if (doMaterial_)
157 material_ = Get<Teuchos::RCP<MultiVector>>(fineLevel, "Material");
158 Teuchos::RCP<LWGraph> fineGraph = Teuchos::null;
159 Teuchos::RCP<LWGraph> coarseGraph = Teuchos::null;
160 if (doFineGraphEdges_)
161 fineGraph = Get<RCP<LWGraph>>(fineLevel, "Graph");
162 if (doCoarseGraphEdges_)
163 coarseGraph = Get<RCP<LWGraph>>(coarseLevel, "Graph");
164 if (useVTK) // otherwise leave null, will not be accessed by non-vtk code
165 {
166 coords = Get<RCP<CoordinateMultiVector>>(fineLevel, "Coordinates");
167 coords_ = coords;
168 if (doCoarseGraphEdges_)
169 coordsCoarse = Get<RCP<CoordinateMultiVector>>(coarseLevel, "Coordinates");
170 dims_ = coords->getNumVectors(); // 2D or 3D?
171 if (numProcs > 1) {
172 if (aggregates->AggregatesCrossProcessors()) { // Do we want to use the map from aggregates here instead of the map from A? Using the map from A seems to be problematic with multiple dofs per node
173 RCP<Import> coordImporter = Xpetra::ImportFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(coords->getMap(), Amat->getColMap());
174 RCP<CoordinateMultiVector> ghostedCoords = Xpetra::MultiVectorFactory<coordinate_type, LocalOrdinal, GlobalOrdinal, Node>::Build(Amat->getColMap(), dims_);
175 ghostedCoords->doImport(*coords, *coordImporter, Xpetra::INSERT);
176 coords = ghostedCoords;
177 coords_ = ghostedCoords;
178 }
179 if (doCoarseGraphEdges_) {
180 RCP<Import> coordImporter = Xpetra::ImportFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(coordsCoarse->getMap(), Ac->getColMap());
181 RCP<CoordinateMultiVector> ghostedCoords = Xpetra::MultiVectorFactory<coordinate_type, LocalOrdinal, GlobalOrdinal, Node>::Build(Ac->getColMap(), dims_);
182 ghostedCoords->doImport(*coordsCoarse, *coordImporter, Xpetra::INSERT);
183 coordsCoarse = ghostedCoords;
184 coordsCoarse_ = ghostedCoords;
185 }
186 }
187 }
188 GetOStream(Runtime0) << "AggregationExportFactory: DofsPerNode: " << DofsPerNode << std::endl;
189
190 Teuchos::RCP<LocalOrdinalMultiVector> vertex2AggId = aggregates->GetVertex2AggId();
191 vertex2AggId_ = vertex2AggId;
192
193 // prepare for calculating global aggregate ids
194 std::vector<GlobalOrdinal> numAggsGlobal(numProcs, 0);
195 std::vector<GlobalOrdinal> numAggsLocal(numProcs, 0);
196 std::vector<GlobalOrdinal> minGlobalAggId(numProcs, 0);
197
198 numAggsLocal[myRank] = aggregates->GetNumAggregates();
199 Teuchos::reduceAll(*comm, Teuchos::REDUCE_SUM, numProcs, &numAggsLocal[0], &numAggsGlobal[0]);
200 for (int i = 1; i < Teuchos::as<int>(numAggsGlobal.size()); ++i) {
201 numAggsGlobal[i] += numAggsGlobal[i - 1];
202 minGlobalAggId[i] = numAggsGlobal[i - 1];
203 }
204 if (numProcs == 0)
205 aggsOffset_ = 0;
206 else
207 aggsOffset_ = minGlobalAggId[myRank];
208 ArrayRCP<LO> aggStart;
209 ArrayRCP<GlobalOrdinal> aggToRowMap;
210 amalgInfo->UnamalgamateAggregates(*aggregates, aggStart, aggToRowMap);
211 int timeStep = pL.get<int>("Output file: time step");
212 int iter = pL.get<int>("Output file: iter");
213 filenameToWrite = this->replaceAll(filenameToWrite, "%LEVELID", toString(fineLevel.GetLevelID()));
214 filenameToWrite = this->replaceAll(filenameToWrite, "%TIMESTEP", toString(timeStep));
215 filenameToWrite = this->replaceAll(filenameToWrite, "%ITER", toString(iter));
216 // Proc id MUST be included in vtu filenames to distinguish them (if multiple procs)
217 // In all other cases (else), including processor # in filename is optional
218 string masterStem = "";
219 if (useVTK) {
220 masterStem = filenameToWrite.substr(0, filenameToWrite.rfind(".vtu"));
221 masterStem = this->replaceAll(masterStem, "%PROCID", "");
222 }
223 pvtuFilename = masterStem + "-master.pvtu";
224 string baseFname = filenameToWrite; // get a version of the filename string with the %PROCID token, but without substituting myRank (needed for pvtu output)
225 filenameToWrite = this->replaceAll(filenameToWrite, "%PROCID", toString(myRank));
226 GetOStream(Runtime0) << "AggregationExportFactory: outputfile \"" << filenameToWrite << "\"" << std::endl;
227 ofstream fout(filenameToWrite.c_str());
228 GO numAggs = aggregates->GetNumAggregates();
229 if (!useVTK) {
230 GO indexBase = aggregates->GetMap()->getIndexBase(); // extract indexBase from overlapping map within aggregates structure. The indexBase is constant throughout the whole simulation (either 0 = C++ or 1 = Fortran)
231 GO offset = amalgInfo->GlobalOffset(); // extract offset for global dof ids
232 vector<GlobalOrdinal> nodeIds;
233 for (int i = 0; i < numAggs; ++i) {
234 fout << "Agg " << minGlobalAggId[myRank] + i << " Proc " << myRank << ":";
235
236 // TODO: Use k+=DofsPerNode instead of ++k and get rid of std::unique call afterwards
237 for (int k = aggStart[i]; k < aggStart[i + 1]; ++k) {
238 nodeIds.push_back((aggToRowMap[k] - offset - indexBase) / DofsPerNode + indexBase);
239 }
240
241 // remove duplicate entries from nodeids
242 std::sort(nodeIds.begin(), nodeIds.end());
243 typename std::vector<GlobalOrdinal>::iterator endLocation = std::unique(nodeIds.begin(), nodeIds.end());
244 nodeIds.erase(endLocation, nodeIds.end());
245
246 // print out nodeids
247 for (typename std::vector<GlobalOrdinal>::iterator printIt = nodeIds.begin(); printIt != nodeIds.end(); printIt++)
248 fout << " " << *printIt;
249 nodeIds.clear();
250 fout << std::endl;
251 }
252 } else {
253 using namespace std;
254 // Make sure we have coordinates
255 TEUCHOS_TEST_FOR_EXCEPTION(coords.is_null(), Exceptions::RuntimeError, "AggExportFactory could not get coordinates, but they are required for VTK output.");
256 numAggs_ = numAggs;
257 numNodes_ = coords->getLocalLength();
258 // Get the sizes of the aggregates to speed up grabbing node IDs
259 aggSizes_ = aggregates->ComputeAggregateSizesArrayRCP();
260 myRank_ = myRank;
261 string aggStyle = "Point Cloud";
262 try {
263 aggStyle = pL.get<string>("aggregation: output file: agg style"); // Let "Point Cloud" be the default style
264 } catch (std::exception& e) {
265 }
266 vector<int> vertices;
267 vector<int> geomSizes;
268 string indent = "";
269 nodeMap_ = Amat->getMap();
270 for (LocalOrdinal i = 0; i < numNodes_; i++) {
271 isRoot_.push_back(aggregates->IsRoot(i));
272 }
273 // If problem is serial and not outputting fine nor coarse graph edges, don't make pvtu file
274 // Otherwise create it
275 if (myRank == 0 && (numProcs != 1 || doCoarseGraphEdges_ || doFineGraphEdges_)) {
276 ofstream pvtu(pvtuFilename.c_str());
277 writePVTU_(pvtu, baseFname, numProcs);
278 pvtu.close();
279 }
280 if (aggStyle == "Point Cloud")
281 this->doPointCloud(vertices, geomSizes, numAggs_, numNodes_);
282 else if (aggStyle == "Jacks") {
283 auto vertex2AggIds = vertex2AggId_->getDataNonConst(0);
284 this->doJacks(vertices, geomSizes, numAggs_, numNodes_, isRoot_, vertex2AggIds);
285 } else if (aggStyle == "Jacks++") // Not actually implemented
286 doJacksPlus_(vertices, geomSizes);
287 else if (aggStyle == "Convex Hulls")
288 doConvexHulls(vertices, geomSizes);
289 else {
290 GetOStream(Warnings0) << " Unrecognized agg style.\n Possible values are Point Cloud, Jacks, Jacks++, and Convex Hulls.\n Defaulting to Point Cloud." << std::endl;
291 aggStyle = "Point Cloud";
292 this->doPointCloud(vertices, geomSizes, numAggs_, numNodes_);
293 }
294 writeFile_(fout, aggStyle, vertices, geomSizes);
295 if (doCoarseGraphEdges_) {
296 string fname = filenameToWrite;
297 string cEdgeFile = fname.insert(fname.rfind(".vtu"), "-coarsegraph");
298 std::ofstream edgeStream(cEdgeFile.c_str());
299 doGraphEdges_(edgeStream, Ac, coarseGraph, false, DofsPerNode);
300 }
301 if (doFineGraphEdges_) {
302 string fname = filenameToWrite;
303 string fEdgeFile = fname.insert(fname.rfind(".vtu"), "-finegraph");
304 std::ofstream edgeStream(fEdgeFile.c_str());
305 doGraphEdges_(edgeStream, Amat, fineGraph, true, DofsPerNode);
306 }
307 if (myRank == 0 && pL.get<bool>("aggregation: output file: build colormap")) {
308 buildColormap_();
309 }
310 }
311 fout.close();
312}
313
314template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
315void AggregationExportFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::doJacksPlus_(std::vector<int>& /* vertices */, std::vector<int>& /* geomSizes */) const {
316 // TODO
317}
318
319template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
320void AggregationExportFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::doConvexHulls(std::vector<int>& vertices, std::vector<int>& geomSizes) const {
321 Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> xCoords = coords_->getData(0);
322 Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> yCoords = coords_->getData(1);
323 Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> zCoords = Teuchos::null;
324
325 auto vertex2AggIds = vertex2AggId_->getDataNonConst(0);
326
327 if (dims_ == 2) {
328 this->doConvexHulls2D(vertices, geomSizes, numAggs_, numNodes_, isRoot_, vertex2AggIds, xCoords, yCoords);
329 } else {
330 zCoords = coords_->getData(2);
331 this->doConvexHulls3D(vertices, geomSizes, numAggs_, numNodes_, isRoot_, vertex2AggIds, xCoords, yCoords, zCoords);
332 }
333}
334
335template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
336void AggregationExportFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::doGraphEdges_(std::ofstream& fout, Teuchos::RCP<Matrix>& A, Teuchos::RCP<LWGraph>& G, bool fine, int dofs) const {
337 using namespace std;
338 ArrayView<const Scalar> values;
339 // Allow two different colors of connections (by setting "aggregates" scalar to CONTRAST_1 or CONTRAST_2)
340 vector<pair<int, int>> vert1; // vertices (node indices)
341 vector<pair<int, int>> vert2; // size of every cell is assumed to be 2 vertices, since all edges are drawn as lines
342
343 Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> xCoords = coords_->getData(0);
344 Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> yCoords = coords_->getData(1);
345 Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> zCoords = Teuchos::null;
346 if (dims_ == 3)
347 zCoords = coords_->getData(2);
348
349 Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> cx = coordsCoarse_->getData(0);
350 Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> cy = coordsCoarse_->getData(1);
351 Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> cz = Teuchos::null;
352 if (dims_ == 3)
353 cz = coordsCoarse_->getData(2);
354
355 if (A->isGloballyIndexed()) {
356 ArrayView<const GlobalOrdinal> indices;
357 for (GlobalOrdinal globRow = 0; globRow < GlobalOrdinal(A->getGlobalNumRows()); globRow++) {
358 if (dofs == 1)
359 A->getGlobalRowView(globRow, indices, values);
360 auto neighbors = G->getNeighborVertices((LocalOrdinal)globRow);
361 int gEdge = 0;
362 int aEdge = 0;
363 while (gEdge != int(neighbors.length)) {
364 if (dofs == 1) {
365 if (neighbors(gEdge) == indices[aEdge]) {
366 // graph and matrix both have this edge, wasn't filtered, show as color 1
367 vert1.push_back(pair<int, int>(int(globRow), neighbors(gEdge)));
368 gEdge++;
369 aEdge++;
370 } else {
371 // graph contains an edge at gEdge which was filtered from A, show as color 2
372 vert2.push_back(pair<int, int>(int(globRow), neighbors(gEdge)));
373 gEdge++;
374 }
375 } else // for multiple DOF problems, don't try to detect filtered edges and ignore A
376 {
377 vert1.push_back(pair<int, int>(int(globRow), neighbors(gEdge)));
378 gEdge++;
379 }
380 }
381 }
382 } else {
383 ArrayView<const LocalOrdinal> indices;
384 for (LocalOrdinal locRow = 0; locRow < LocalOrdinal(A->getLocalNumRows()); locRow++) {
385 if (dofs == 1)
386 A->getLocalRowView(locRow, indices, values);
387 auto neighbors = G->getNeighborVertices(locRow);
388 // Add those local indices (columns) to the list of connections (which are pairs of the form (localM, localN))
389 int gEdge = 0;
390 int aEdge = 0;
391 while (gEdge != int(neighbors.length)) {
392 if (dofs == 1) {
393 if (neighbors(gEdge) == indices[aEdge]) {
394 vert1.push_back(pair<int, int>(locRow, neighbors(gEdge)));
395 gEdge++;
396 aEdge++;
397 } else {
398 vert2.push_back(pair<int, int>(locRow, neighbors(gEdge)));
399 gEdge++;
400 }
401 } else {
402 vert1.push_back(pair<int, int>(locRow, neighbors(gEdge)));
403 gEdge++;
404 }
405 }
406 }
407 }
408 for (size_t i = 0; i < vert1.size(); i++) {
409 if (vert1[i].first > vert1[i].second) {
410 int temp = vert1[i].first;
411 vert1[i].first = vert1[i].second;
412 vert1[i].second = temp;
413 }
414 }
415 for (size_t i = 0; i < vert2.size(); i++) {
416 if (vert2[i].first > vert2[i].second) {
417 int temp = vert2[i].first;
418 vert2[i].first = vert2[i].second;
419 vert2[i].second = temp;
420 }
421 }
422 sort(vert1.begin(), vert1.end());
423 vector<pair<int, int>>::iterator newEnd = unique(vert1.begin(), vert1.end()); // remove duplicate edges
424 vert1.erase(newEnd, vert1.end());
425 sort(vert2.begin(), vert2.end());
426 newEnd = unique(vert2.begin(), vert2.end());
427 vert2.erase(newEnd, vert2.end());
428 vector<int> points1;
429 points1.reserve(2 * vert1.size());
430 for (size_t i = 0; i < vert1.size(); i++) {
431 points1.push_back(vert1[i].first);
432 points1.push_back(vert1[i].second);
433 }
434 vector<int> points2;
435 points2.reserve(2 * vert2.size());
436 for (size_t i = 0; i < vert2.size(); i++) {
437 points2.push_back(vert2[i].first);
438 points2.push_back(vert2[i].second);
439 }
440 vector<int> unique1 = this->makeUnique(points1);
441 vector<int> unique2 = this->makeUnique(points2);
442 fout << "<VTKFile type=\"UnstructuredGrid\" byte_order=\"LittleEndian\">" << endl;
443 fout << " <UnstructuredGrid>" << endl;
444 fout << " <Piece NumberOfPoints=\"" << unique1.size() + unique2.size() << "\" NumberOfCells=\"" << vert1.size() + vert2.size() << "\">" << endl;
445 fout << " <PointData Scalars=\"Node Aggregate Processor\">" << endl;
446 fout << " <DataArray type=\"Int32\" Name=\"Node\" format=\"ascii\">" << endl; // node and aggregate will be set to CONTRAST_1|2, but processor will have its actual value
447 string indent = " ";
448 fout << indent;
449 for (size_t i = 0; i < unique1.size(); i++) {
450 fout << CONTRAST_1_ << " ";
451 if (i % 25 == 24)
452 fout << endl
453 << indent;
454 }
455 for (size_t i = 0; i < unique2.size(); i++) {
456 fout << CONTRAST_2_ << " ";
457 if ((i + 2 * vert1.size()) % 25 == 24)
458 fout << endl
459 << indent;
460 }
461 fout << endl;
462 fout << " </DataArray>" << endl;
463 fout << " <DataArray type=\"Int32\" Name=\"Aggregate\" format=\"ascii\">" << endl;
464 fout << indent;
465 for (size_t i = 0; i < unique1.size(); i++) {
466 fout << CONTRAST_1_ << " ";
467 if (i % 25 == 24)
468 fout << endl
469 << indent;
470 }
471 for (size_t i = 0; i < unique2.size(); i++) {
472 fout << CONTRAST_2_ << " ";
473 if ((i + 2 * vert2.size()) % 25 == 24)
474 fout << endl
475 << indent;
476 }
477 fout << endl;
478 fout << " </DataArray>" << endl;
479 fout << " <DataArray type=\"Int32\" Name=\"Processor\" format=\"ascii\">" << endl;
480 fout << indent;
481 for (size_t i = 0; i < unique1.size() + unique2.size(); i++) {
482 fout << myRank_ << " ";
483 if (i % 25 == 24)
484 fout << endl
485 << indent;
486 }
487 fout << endl;
488 fout << " </DataArray>" << endl;
489 fout << " </PointData>" << endl;
490 fout << " <Points>" << endl;
491 fout << " <DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"ascii\">" << endl;
492 fout << indent;
493 for (size_t i = 0; i < unique1.size(); i++) {
494 if (fine) {
495 fout << xCoords[unique1[i]] << " " << yCoords[unique1[i]] << " ";
496 if (dims_ == 3)
497 fout << zCoords[unique1[i]] << " ";
498 else
499 fout << "0 ";
500 if (i % 2)
501 fout << endl
502 << indent;
503 } else {
504 fout << cx[unique1[i]] << " " << cy[unique1[i]] << " ";
505 if (dims_ == 3)
506 fout << cz[unique1[i]] << " ";
507 else
508 fout << "0 ";
509 if (i % 2)
510 fout << endl
511 << indent;
512 }
513 }
514 for (size_t i = 0; i < unique2.size(); i++) {
515 if (fine) {
516 fout << xCoords[unique2[i]] << " " << yCoords[unique2[i]] << " ";
517 if (dims_ == 3)
518 fout << zCoords[unique2[i]] << " ";
519 else
520 fout << "0 ";
521 if (i % 2)
522 fout << endl
523 << indent;
524 } else {
525 fout << cx[unique2[i]] << " " << cy[unique2[i]] << " ";
526 if (dims_ == 3)
527 fout << cz[unique2[i]] << " ";
528 else
529 fout << "0 ";
530 if ((i + unique1.size()) % 2)
531 fout << endl
532 << indent;
533 }
534 }
535 fout << endl;
536 fout << " </DataArray>" << endl;
537 fout << " </Points>" << endl;
538 fout << " <Cells>" << endl;
539 fout << " <DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">" << endl;
540 fout << indent;
541 for (size_t i = 0; i < points1.size(); i++) {
542 fout << points1[i] << " ";
543 if (i % 10 == 9)
544 fout << endl
545 << indent;
546 }
547 for (size_t i = 0; i < points2.size(); i++) {
548 fout << points2[i] + unique1.size() << " ";
549 if ((i + points1.size()) % 10 == 9)
550 fout << endl
551 << indent;
552 }
553 fout << endl;
554 fout << " </DataArray>" << endl;
555 fout << " <DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">" << endl;
556 fout << indent;
557 int offset = 0;
558 for (size_t i = 0; i < vert1.size() + vert2.size(); i++) {
559 offset += 2;
560 fout << offset << " ";
561 if (i % 10 == 9)
562 fout << endl
563 << indent;
564 }
565 fout << endl;
566 fout << " </DataArray>" << endl;
567 fout << " <DataArray type=\"Int32\" Name=\"types\" format=\"ascii\">" << endl;
568 fout << indent;
569 for (size_t i = 0; i < vert1.size() + vert2.size(); i++) {
570 fout << "3 ";
571 if (i % 25 == 24)
572 fout << endl
573 << indent;
574 }
575 fout << endl;
576 fout << " </DataArray>" << endl;
577 fout << " </Cells>" << endl;
578 fout << " </Piece>" << endl;
579 fout << " </UnstructuredGrid>" << endl;
580 fout << "</VTKFile>" << endl;
581}
582
583template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
584void AggregationExportFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeFile_(std::ofstream& fout, std::string styleName, std::vector<int>& vertices, std::vector<int>& geomSizes) const {
585 using namespace std;
586
587 Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> xCoords = coords_->getData(0);
588 Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> yCoords = coords_->getData(1);
589 Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> zCoords = Teuchos::null;
590 if (dims_ == 3)
591 zCoords = coords_->getData(2);
592
593 auto vertex2AggIds = vertex2AggId_->getDataNonConst(0);
594
595 Teuchos::ArrayRCP<const Scalar> qualities;
596 if (doAggQuality_)
597 qualities = qualities_->getData(0);
598
599 Teuchos::ArrayRCP<Teuchos::ArrayRCP<const Scalar>> material;
600 if (doMaterial_) {
601 size_t dim = material_->getNumVectors();
602 material.resize(dim);
603 for (size_t k = 0; k < dim; k++)
604 material[k] = material_->getData(k);
605 }
606
607 vector<int> uniqueFine = this->makeUnique(vertices);
608 string indent = " ";
609 fout << "<!--" << styleName << " Aggregates Visualization-->" << endl;
610 fout << "<VTKFile type=\"UnstructuredGrid\" byte_order=\"LittleEndian\">" << endl;
611 fout << " <UnstructuredGrid>" << endl;
612 fout << " <Piece NumberOfPoints=\"" << uniqueFine.size() << "\" NumberOfCells=\"" << geomSizes.size() << "\">" << endl;
613 fout << " <PointData Scalars=\"Node Aggregate Processor\">" << endl;
614 fout << " <DataArray type=\"Int32\" Name=\"Node\" format=\"ascii\">" << endl;
615 indent = " ";
616 fout << indent;
617 bool localIsGlobal = GlobalOrdinal(nodeMap_->getGlobalNumElements()) == GlobalOrdinal(nodeMap_->getLocalNumElements());
618 for (size_t i = 0; i < uniqueFine.size(); i++) {
619 if (localIsGlobal) {
620 fout << uniqueFine[i] << " "; // if all nodes are on this processor, do not map from local to global
621 } else
622 fout << nodeMap_->getGlobalElement(uniqueFine[i]) << " ";
623 if (i % 10 == 9)
624 fout << endl
625 << indent;
626 }
627 fout << endl;
628 fout << " </DataArray>" << endl;
629 fout << " <DataArray type=\"Int32\" Name=\"Aggregate\" format=\"ascii\">" << endl;
630 fout << indent;
631 for (size_t i = 0; i < uniqueFine.size(); i++) {
632 if (vertex2AggIds[uniqueFine[i]] == -1)
633 fout << vertex2AggIds[uniqueFine[i]] << " ";
634 else
635 fout << aggsOffset_ + vertex2AggIds[uniqueFine[i]] << " ";
636 if (i % 10 == 9)
637 fout << endl
638 << indent;
639 }
640 fout << endl;
641 fout << " </DataArray>" << endl;
642 fout << " <DataArray type=\"Int32\" Name=\"Processor\" format=\"ascii\">" << endl;
643 fout << indent;
644 for (size_t i = 0; i < uniqueFine.size(); i++) {
645 fout << myRank_ << " ";
646 if (i % 20 == 19)
647 fout << endl
648 << indent;
649 }
650 fout << endl;
651 fout << " </DataArray>" << endl;
652 if (doAggQuality_) {
653 fout << " <DataArray type=\"Float64\" Name=\"Quality\" format=\"ascii\">" << endl;
654 fout << indent;
655 for (size_t i = 0; i < uniqueFine.size(); i++) {
656 fout << qualities[vertex2AggIds[uniqueFine[i]]] << " ";
657 if (i % 10 == 9)
658 fout << endl
659 << indent;
660 }
661 fout << endl;
662 fout << " </DataArray>" << endl;
663 }
664 // Material stuff
665 if (doMaterial_) {
666 size_t dim = material_->getNumVectors();
667 fout << " <DataArray type=\"Float64\" NumberOfComponents=\"" << dim << "\" Name=\"Material\" format=\"ascii\">" << endl;
668 fout << indent;
669 for (size_t i = 0; i < uniqueFine.size(); i++) {
670 for (size_t k = 0; k < dim; k++) {
671 fout << material[k][uniqueFine[i]] << " ";
672 }
673 fout << endl
674 << indent;
675 }
676 fout << endl;
677 fout << " </DataArray>" << endl;
678 }
679 fout << " </PointData>" << endl;
680 fout << " <Points>" << endl;
681 fout << " <DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"ascii\">" << endl;
682 fout << indent;
683 for (size_t i = 0; i < uniqueFine.size(); i++) {
684 fout << xCoords[uniqueFine[i]] << " " << yCoords[uniqueFine[i]] << " ";
685 if (dims_ == 2)
686 fout << "0 ";
687 else
688 fout << zCoords[uniqueFine[i]] << " ";
689 if (i % 3 == 2)
690 fout << endl
691 << indent;
692 }
693 fout << endl;
694 fout << " </DataArray>" << endl;
695 fout << " </Points>" << endl;
696 fout << " <Cells>" << endl;
697 fout << " <DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">" << endl;
698 fout << indent;
699 for (size_t i = 0; i < vertices.size(); i++) {
700 fout << vertices[i] << " ";
701 if (i % 10 == 9)
702 fout << endl
703 << indent;
704 }
705 fout << endl;
706 fout << " </DataArray>" << endl;
707 fout << " <DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">" << endl;
708 fout << indent;
709 int accum = 0;
710 for (size_t i = 0; i < geomSizes.size(); i++) {
711 accum += geomSizes[i];
712 fout << accum << " ";
713 if (i % 10 == 9)
714 fout << endl
715 << indent;
716 }
717 fout << endl;
718 fout << " </DataArray>" << endl;
719 fout << " <DataArray type=\"Int32\" Name=\"types\" format=\"ascii\">" << endl;
720 fout << indent;
721 for (size_t i = 0; i < geomSizes.size(); i++) {
722 switch (geomSizes[i]) {
723 case 1:
724 fout << "1 "; // Point
725 break;
726 case 2:
727 fout << "3 "; // Line
728 break;
729 case 3:
730 fout << "5 "; // Triangle
731 break;
732 default:
733 fout << "7 "; // Polygon
734 }
735 if (i % 30 == 29)
736 fout << endl
737 << indent;
738 }
739 fout << endl;
740 fout << " </DataArray>" << endl;
741 fout << " </Cells>" << endl;
742 fout << " </Piece>" << endl;
743 fout << " </UnstructuredGrid>" << endl;
744 fout << "</VTKFile>" << endl;
745}
746
747template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
748void AggregationExportFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::buildColormap_() const {
749 using namespace std;
750 try {
751 ofstream color("random-colormap.xml");
752 color << "<ColorMap name=\"MueLu-Random\" space=\"RGB\">" << endl;
753 // Give -1, -2, -3 distinctive colors (so that the style functions can have constrasted geometry types)
754 // Do red, orange, yellow to constrast with cool color spectrum for other types
755 color << " <Point x=\"" << CONTRAST_1_ << "\" o=\"1\" r=\"1\" g=\"0\" b=\"0\"/>" << endl;
756 color << " <Point x=\"" << CONTRAST_2_ << "\" o=\"1\" r=\"1\" g=\"0.6\" b=\"0\"/>" << endl;
757 color << " <Point x=\"" << CONTRAST_3_ << "\" o=\"1\" r=\"1\" g=\"1\" b=\"0\"/>" << endl;
758 srand(time(NULL));
759 for (int i = 0; i < 5000; i += 4) {
760 color << " <Point x=\"" << i << "\" o=\"1\" r=\"" << (rand() % 50) / 256.0 << "\" g=\"" << (rand() % 256) / 256.0 << "\" b=\"" << (rand() % 256) / 256.0 << "\"/>" << endl;
761 }
762 color << "</ColorMap>" << endl;
763 color.close();
764 } catch (std::exception& e) {
765 GetOStream(Warnings0) << " Error while building colormap file: " << e.what() << endl;
766 }
767}
768
769template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
770void AggregationExportFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writePVTU_(std::ofstream& pvtu, std::string baseFname, int numProcs) const {
771 using namespace std;
772 // If using vtk, filenameToWrite now contains final, correct ***.vtu filename (for the current proc)
773 // So the root proc will need to use its own filenameToWrite to make a list of the filenames of all other procs to put in
774 // pvtu file.
775 pvtu << "<VTKFile type=\"PUnstructuredGrid\" byte_order=\"LittleEndian\">" << endl;
776 pvtu << " <PUnstructuredGrid GhostLevel=\"0\">" << endl;
777 pvtu << " <PPointData Scalars=\"Node Aggregate Processor\">" << endl;
778 pvtu << " <PDataArray type=\"Int32\" Name=\"Node\"/>" << endl;
779 pvtu << " <PDataArray type=\"Int32\" Name=\"Aggregate\"/>" << endl;
780 pvtu << " <PDataArray type=\"Int32\" Name=\"Processor\"/>" << endl;
781 pvtu << " </PPointData>" << endl;
782 pvtu << " <PPoints>" << endl;
783 pvtu << " <PDataArray type=\"Float64\" NumberOfComponents=\"3\"/>" << endl;
784 pvtu << " </PPoints>" << endl;
785 for (int i = 0; i < numProcs; i++) {
786 // specify the piece for each proc (the replaceAll expression matches what the filenames will be on other procs)
787 pvtu << " <Piece Source=\"" << this->replaceAll(baseFname, "%PROCID", toString(i)) << "\"/>" << endl;
788 }
789 // reference files with graph pieces, if applicable
790 if (doFineGraphEdges_) {
791 for (int i = 0; i < numProcs; i++) {
792 string fn = this->replaceAll(baseFname, "%PROCID", toString(i));
793 pvtu << " <Piece Source=\"" << fn.insert(fn.rfind(".vtu"), "-finegraph") << "\"/>" << endl;
794 }
795 }
796 if (doCoarseGraphEdges_) {
797 for (int i = 0; i < numProcs; i++) {
798 string fn = this->replaceAll(baseFname, "%PROCID", toString(i));
799 pvtu << " <Piece Source=\"" << fn.insert(fn.rfind(".vtu"), "-coarsegraph") << "\"/>" << endl;
800 }
801 }
802 pvtu << " </PUnstructuredGrid>" << endl;
803 pvtu << "</VTKFile>" << endl;
804 pvtu.close();
805}
806} // namespace MueLu
807
808#endif /* MUELU_AGGREGATIONEXPORTFACTORY_DEF_HPP_ */
MueLu_AggregationExportFactory_decl.hpp
MueLu_FactoryManagerBase.hpp
MueLu_Level.hpp
MueLu_Monitor.hpp
LocalOrdinal
MueLu::DefaultLocalOrdinal LocalOrdinal
Definition MueLu_UseDefaultTypes.hpp:13
GlobalOrdinal
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Definition MueLu_UseDefaultTypes.hpp:14
MueLu::AggregationExportFactory::doJacksPlus_
void doJacksPlus_(std::vector< int > &vertices, std::vector< int > &geomSizes) const
Definition MueLu_AggregationExportFactory_def.hpp:315
MueLu::AggregationExportFactory::writePVTU_
void writePVTU_(std::ofstream &pvtu, std::string baseFname, int numProcs) const
Definition MueLu_AggregationExportFactory_def.hpp:770
MueLu::AggregationExportFactory::doConvexHulls
void doConvexHulls(std::vector< int > &vertices, std::vector< int > &geomSizes) const
Definition MueLu_AggregationExportFactory_def.hpp:320
MueLu::AggregationExportFactory::~AggregationExportFactory
virtual ~AggregationExportFactory()
Destructor.
MueLu::AggregationExportFactory::doGraphEdges_
void doGraphEdges_(std::ofstream &fout, Teuchos::RCP< Matrix > &A, Teuchos::RCP< LWGraph > &G, bool fine, int dofs) const
Definition MueLu_AggregationExportFactory_def.hpp:336
MueLu::AggregationExportFactory::writeFile_
void writeFile_(std::ofstream &fout, std::string styleName, std::vector< int > &vertices, std::vector< int > &geomSizes) const
Definition MueLu_AggregationExportFactory_def.hpp:584
MueLu::AggregationExportFactory::Build
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
Definition MueLu_AggregationExportFactory_def.hpp:119
MueLu::AggregationExportFactory::GetValidParameterList
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Definition MueLu_AggregationExportFactory_def.hpp:56
MueLu::AggregationExportFactory::DeclareInput
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
Definition MueLu_AggregationExportFactory_def.hpp:91
MueLu::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition MueLu_Exceptions.hpp:36
MueLu::FactoryMonitor
Timer to be used in factories. Similar to Monitor but with additional timers.
Definition MueLu_Monitor.hpp:165
MueLu::Level
Class that holds all level-specific information.
Definition MueLu_Level.hpp:63
MueLu::Level::GetLevelID
int GetLevelID() const
Return level number.
Definition MueLu_Level.cpp:51
MueLu
Namespace for MueLu classes and methods.
Definition MueLu_BrickAggregationFactory_decl.hpp:42
MueLu::Warnings0
@ Warnings0
Important warning messages (one line)
Definition MueLu_VerbosityLevel.hpp:20
MueLu::Runtime0
@ Runtime0
One-liner description of what is happening.
Definition MueLu_VerbosityLevel.hpp:25
MueLu::replaceAll
void replaceAll(std::string &str, const std::string &from, const std::string &to)
Definition MueLu_ParameterListUtils.cpp:134
MueLu::toString
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
Definition MueLu_Utilities_decl.hpp:599