10 #ifndef MUELU_COARSENINGVISUALIZATIONFACTORY_DEF_HPP_
11 #define MUELU_COARSENINGVISUALIZATIONFACTORY_DEF_HPP_
15 #include <Xpetra_MultiVectorFactory.hpp>
21 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
25 validParamList->
set<
int>(
"visualization: start level", 0,
"visualize only levels with level ids greater or equal than start level");
27 validParamList->
set<
RCP<const FactoryBase> >(
"P", Teuchos::null,
"Prolongator factory. The user has to declare either P or Ptent but not both at the same time.");
28 validParamList->
set<
RCP<const FactoryBase> >(
"Ptent", Teuchos::null,
"Tentative prolongator factory. The user has to declare either P or Ptent as input but not both at the same time");
32 return validParamList;
35 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
37 this->Input(fineLevel,
"Coordinates");
41 pL.
isParameter(
"Ptent") && GetFactory(
"Ptent") != Teuchos::null,
43 "You must not declare both P and Ptent. Use only once for visualization.");
45 "You have to either declare P or Ptent for visualization, but not both.");
47 if (GetFactory(
"P") != Teuchos::null && GetFactory(
"Ptent") == Teuchos::null)
48 this->Input(coarseLevel,
"P");
49 else if (GetFactory(
"Ptent") != Teuchos::null && GetFactory(
"P") == Teuchos::null)
50 this->Input(coarseLevel,
"Ptent");
52 if (pL.
get<
bool>(
"visualization: fine graph edges"))
53 Input(fineLevel,
"Graph");
55 if(pL.
get<
bool>(
"visualization: coarse graph edges")) {
56 Input(coarseLevel,
"Coordinates");
57 Input(coarseLevel,
"Graph");
62 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
67 if (this->GetFactory(
"P") != Teuchos::null && this->GetFactory(
"Ptent") == Teuchos::null)
68 P = Get<RCP<Matrix> >(coarseLevel,
"P");
69 if (GetFactory(
"Ptent") != Teuchos::null && GetFactory(
"P") == Teuchos::null)
77 if (P->IsView(
"stridedMaps") && Teuchos::rcp_dynamic_cast<
const StridedMap>(P->getRowMap(
"stridedMaps")) != Teuchos::null) {
78 strRowMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(P->getRowMap(
"stridedMaps"));
81 std::vector<size_t> stridingInfo = strRowMap->getStridingData();
82 for (
size_t j = 0; j < Teuchos::as<size_t>(blockid); j++)
83 stridedRowOffset += stridingInfo[j];
84 dofsPerNode = Teuchos::as<LocalOrdinal>(stridingInfo[blockid]);
86 dofsPerNode = strRowMap->getFixedBlockSize();
88 GetOStream(
Runtime1) <<
"CoarseningVisualizationFactory::Build():"
89 <<
" #dofs per node = " << dofsPerNode << std::endl;
95 if (P->IsView(
"stridedMaps") && Teuchos::rcp_dynamic_cast<
const StridedMap>(P->getRowMap(
"stridedMaps")) != Teuchos::null) {
96 strDomainMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(P->getColMap(
"stridedMaps"));
97 LocalOrdinal blockid = strDomainMap->getStridedBlockId();
100 std::vector<size_t> stridingInfo = strDomainMap->getStridingData();
101 for (
size_t j = 0; j < Teuchos::as<size_t>(blockid); j++)
102 stridedColumnOffset += stridingInfo[j];
103 columnsPerNode = Teuchos::as<LocalOrdinal>(stridingInfo[blockid]);
105 columnsPerNode = strDomainMap->getFixedBlockSize();
107 GetOStream(
Runtime1) <<
"CoarseningVisualizationFactory::Build():"
108 <<
" #columns per node = " << columnsPerNode << std::endl;
113 "CoarseningVisualization only supports non-overlapping transfers");
116 LocalOrdinal numLocalAggs = strDomainMap->getLocalNumElements() / columnsPerNode;
117 std::vector<std::set<LocalOrdinal> > localAggs(numLocalAggs);
120 for (
LO row = 0; row < Teuchos::as<LO>(P->getRowMap()->getLocalNumElements()); ++row) {
123 P->getLocalRowView(row, indices, vals);
126 localAggs[(*c) / columnsPerNode].insert(row / dofsPerNode);
131 std::vector<int> myLocalAggsPerProc(comm->getSize(), 0);
132 std::vector<int> numLocalAggsPerProc(comm->getSize(), 0);
133 myLocalAggsPerProc[comm->getRank()] = numLocalAggs;
134 Teuchos::reduceAll(*comm,
Teuchos::REDUCE_MAX, comm->getSize(), &myLocalAggsPerProc[0], &numLocalAggsPerProc[0]);
137 for (
int i = 0; i < comm->getRank(); ++i) {
138 myAggOffset += numLocalAggsPerProc[i];
154 "Number of fine level nodes in coordinates is inconsistent with dof based information");
157 if (pL.
get<
bool>(
"visualization: fine graph edges")) {
158 fineGraph = Get<RCP<LWGraph> >(fineLevel,
"Graph");
163 coords = ghostedCoords;
171 if (coords->getNumVectors() == 3) {
176 LocalOrdinal numFineNodes = Teuchos::as<LocalOrdinal>(coords->getLocalLength());
180 for (LocalOrdinal i = 0; i < numLocalAggs; ++i) {
182 for (
typename std::set<LocalOrdinal>::iterator it = localAggs[i].begin(); it != localAggs[i].end(); ++it) {
183 vertex2AggId[*it] = i;
190 std::vector<bool> isRoot(numFineNodes,
false);
191 for (LocalOrdinal i = 0; i < numLocalAggs; ++i) {
197 for (
typename std::set<LocalOrdinal>::iterator it = localAggs[i].begin(); it != localAggs[i].end(); ++it) {
198 xCenter += xCoords[*it];
199 yCenter += yCoords[*it];
200 if (coords->getNumVectors() == 3) zCenter += zCoords[*it];
202 xCenter /= Teuchos::as<LocalOrdinal>(localAggs[i].size());
203 yCenter /= Teuchos::as<LocalOrdinal>(localAggs[i].size());
204 zCenter /= Teuchos::as<LocalOrdinal>(localAggs[i].size());
207 LocalOrdinal rootCandidate = -1;
209 for (
typename std::set<LocalOrdinal>::iterator it = localAggs[i].begin(); it != localAggs[i].end(); ++it) {
213 if (coords->getNumVectors() == 3) tempz = zCenter - zCoords[*it];
215 mydistance += tempx * tempx;
216 mydistance += tempy * tempy;
217 mydistance += tempz * tempz;
218 mydistance = sqrt(mydistance);
219 if (mydistance <= minDistance) {
220 minDistance = mydistance;
225 isRoot[rootCandidate] =
true;
228 std::vector<LocalOrdinal> vertices;
229 std::vector<LocalOrdinal> geomSize;
230 int vizLevel = pL.
get<
int>(
"visualization: start level");
231 if (vizLevel <= fineLevel.GetLevelID()) {
232 std::string aggStyle = pL.
get<std::string>(
"visualization: style");
233 if (aggStyle ==
"Point Cloud")
234 this->doPointCloud(vertices, geomSize, numLocalAggs, numFineNodes);
235 else if (aggStyle ==
"Jacks")
236 this->doJacks(vertices, geomSize, numLocalAggs, numFineNodes, isRoot, vertex2AggId);
237 else if (aggStyle ==
"Convex Hulls") {
240 if (coords->getNumVectors() == 3)
241 this->doConvexHulls3D(vertices, geomSize, numLocalAggs, numFineNodes, isRoot, vertex2AggId, xCoords, yCoords, zCoords);
242 else if (coords->getNumVectors() == 2)
243 this->doConvexHulls2D(vertices, geomSize, numLocalAggs, numFineNodes, isRoot, vertex2AggId, xCoords, yCoords);
245 GetOStream(
Warnings0) <<
" Warning: Unrecognized agg style.\nPossible values are Point Cloud, Jacks, Convex Hulls.\nDefaulting to Point Cloud." << std::endl;
246 aggStyle =
"Point Cloud";
247 this->doPointCloud(vertices, geomSize, numLocalAggs, numFineNodes);
252 if (pL.
get<
bool>(
"visualization: fine graph edges")) {
254 "Could not get information about fine graph.");
256 std::vector<LocalOrdinal> fine_edges_vertices;
257 std::vector<LocalOrdinal> fine_edges_geomSize;
258 this->doGraphEdges(fine_edges_vertices, fine_edges_geomSize, fineGraph, xCoords, yCoords, zCoords);
260 std::string fEdgeFineFile = this->getFileName(comm->getSize(), comm->getRank(), fineLevel.GetLevelID(), pL);
261 std::string fEdgeFile = fEdgeFineFile.insert(fEdgeFineFile.rfind(
".vtu"),
"-finegraph");
262 std::ofstream edgeStream(fEdgeFile.c_str());
264 std::vector<int> uniqueFineEdges = this->makeUnique(fine_edges_vertices);
265 this->writeFileVTKOpening(edgeStream, uniqueFineEdges, fine_edges_geomSize);
266 this->writeFileVTKNodes(edgeStream, uniqueFineEdges, nodeMap);
267 this->writeFileVTKData(edgeStream, uniqueFineEdges, myAggOffset, vertex2AggId, comm->getRank());
268 this->writeFileVTKCoordinates(edgeStream, uniqueFineEdges, xCoords, yCoords, zCoords, coords->getNumVectors());
269 this->writeFileVTKCells(edgeStream, uniqueFineEdges, fine_edges_vertices, fine_edges_geomSize);
270 this->writeFileVTKClosing(edgeStream);
275 #if 0 // we don't have access to the coarse graph
276 if (pL.
get<
bool>(
"visualization: coarse graph edges")) {
277 RCP<LWGraph> coarseGraph = Get<RCP<LWGraph> >(coarseLevel,
"Graph");
283 coarseghostedCoords->doImport(*coarsecoords, *coarsecoordImporter,
Xpetra::INSERT);
284 coarsecoords = coarseghostedCoords;
289 if(coarsecoords->getNumVectors() == 3) {
295 std::vector<LocalOrdinal> coarse_edges_vertices;
296 std::vector<LocalOrdinal> coarse_edges_geomSize;
297 this->doGraphEdges(coarse_edges_vertices, coarse_edges_geomSize, coarseGraph, cx, cy, cz);
299 std::string cEdgeFineFile = this->getFileName(comm->getSize(), comm->getRank(), coarseLevel.GetLevelID(), pL);
300 std::string cEdgeFile = cEdgeFineFile.insert(cEdgeFineFile.rfind(
".vtu"),
"-coarsegraph");
301 std::ofstream cedgeStream(cEdgeFile.c_str());
303 std::vector<int> uniqueCoarseEdges = this->makeUnique(coarse_edges_vertices);
304 this->writeFileVTKOpening(cedgeStream, uniqueCoarseEdges, coarse_edges_geomSize);
305 this->writeFileVTKNodes(cedgeStream, uniqueCoarseEdges, coarsenodeMap);
307 this->writeFileVTKCoordinates(cedgeStream, uniqueCoarseEdges, cx, cy, cz, coarsecoords->getNumVectors());
308 this->writeFileVTKCells(cedgeStream, uniqueCoarseEdges, coarse_edges_vertices, coarse_edges_geomSize);
309 this->writeFileVTKClosing(cedgeStream);
314 if (pL.
get<
int>(
"visualization: start level") <= fineLevel.GetLevelID()) {
316 std::string filenameToWrite = this->getFileName(comm->getSize(), comm->getRank(), fineLevel.GetLevelID(), pL);
317 std::ofstream fout(filenameToWrite.c_str());
319 std::vector<int> uniqueFine = this->makeUnique(vertices);
320 this->writeFileVTKOpening(fout, uniqueFine, geomSize);
321 this->writeFileVTKNodes(fout, uniqueFine, nodeMap);
322 this->writeFileVTKData(fout, uniqueFine, myAggOffset, vertex2AggId, comm->getRank());
323 this->writeFileVTKCoordinates(fout, uniqueFine, xCoords, yCoords, zCoords, coords->getNumVectors());
324 this->writeFileVTKCells(fout, uniqueFine, vertices, geomSize);
325 this->writeFileVTKClosing(fout);
329 if (comm->getRank() == 0) {
330 std::string pvtuFilename = this->getPVTUFileName(comm->getSize(), comm->getRank(), fineLevel.GetLevelID(), pL);
331 std::string baseFname = this->getBaseFileName(comm->getSize(), fineLevel.GetLevelID(), pL);
332 std::ofstream pvtu(pvtuFilename.c_str());
333 this->writePVTU(pvtu, baseFname, comm->getSize(), pL.
get<
bool>(
"visualization: fine graph edges"));
338 if (comm->getRank() == 0 && pL.
get<
bool>(
"visualization: build colormap")) {
339 this->buildColormap();
Important warning messages (one line)
MueLu::DefaultLocalOrdinal LocalOrdinal
T & get(const std::string &name, T def_value)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
ParameterList & set(std::string const &name, T &&value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
RCP< ParameterList > GetValidParameterList() const
bool isParameter(const std::string &name) const
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Class that holds all level-specific information.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
const RCP< const Map > GetImportMap() const
Return overlapping import map (nodes).
Exception throws to report errors in the internal logical of the program.
Description of what is happening (more verbose)
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
static RCP< Import< LocalOrdinal, GlobalOrdinal, Node > > Build(const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &source, const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &target, const Teuchos::RCP< Teuchos::ParameterList > &plist=Teuchos::null)