MueLu  Version of the Day
MueLu_VisualizationHelpers_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // MueLu: A package for multigrid based preconditioning
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact
39 // Jonathan Hu (jhu@sandia.gov)
40 // Andrey Prokopenko (aprokop@sandia.gov)
41 // Ray Tuminaro (rstumin@sandia.gov)
42 //
43 // ***********************************************************************
44 //
45 // @HEADER
46 
47 #ifndef MUELU_VISUALIZATIONHELPERS_DEF_HPP_
48 #define MUELU_VISUALIZATIONHELPERS_DEF_HPP_
49 
50 #include <Xpetra_Matrix.hpp>
51 #include <Xpetra_CrsMatrixWrap.hpp>
52 #include <Xpetra_ImportFactory.hpp>
53 #include <Xpetra_MultiVectorFactory.hpp>
55 #include "MueLu_Level.hpp"
56 #include "MueLu_Graph.hpp"
57 #include "MueLu_Monitor.hpp"
58 #include <vector>
59 #include <list>
60 #include <algorithm>
61 #include <string>
62 
63 namespace MueLu {
64 
65  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
67  RCP<ParameterList> validParamList = rcp(new ParameterList());
68 
69  validParamList->set< std::string > ("visualization: output filename", "viz%LEVELID", "filename for VTK-style visualization output");
70  validParamList->set< int > ("visualization: output file: time step", 0, "time step variable for output file name");// Remove me?
71  validParamList->set< int > ("visualization: output file: iter", 0, "nonlinear iteration variable for output file name");//Remove me?
72  validParamList->set<std::string> ("visualization: style", "Point Cloud", "style of aggregate visualization for VTK output. Can be 'Point Cloud', 'Jacks', 'Convex Hulls'");
73  validParamList->set<bool> ("visualization: build colormap", false, "Whether to build a random color map in a separate xml file.");
74  validParamList->set<bool> ("visualization: fine graph edges", false, "Whether to draw all fine node connections along with the aggregates.");
75 
76  return validParamList;
77  }
78 
79  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
80  void VisualizationHelpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::doPointCloud(std::vector<int>& vertices, std::vector<int>& geomSizes, LO numLocalAggs, LO numFineNodes) {
81  vertices.reserve(numFineNodes);
82  geomSizes.reserve(numFineNodes);
83  for(LocalOrdinal i = 0; i < numFineNodes; i++)
84  {
85  vertices.push_back(i);
86  geomSizes.push_back(1);
87  }
88  }
89 
90  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
91  void VisualizationHelpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::doJacks(std::vector<int>& vertices, std::vector<int>& geomSizes, LO numLocalAggs, LO numFineNodes, const std::vector<bool>& isRoot, const ArrayRCP<LO>& vertex2AggId) {
92  //For each aggregate, find the root node then connect it with the other nodes in the aggregate
93  //Doesn't matter the order, as long as all edges are found.
94  vertices.reserve(vertices.size() + 3 * (numFineNodes - numLocalAggs));
95  geomSizes.reserve(vertices.size() + 2 * (numFineNodes - numLocalAggs));
96  int root = 0;
97  for(int i = 0; i < numLocalAggs; i++) //TODO: Replace this O(n^2) with a better way
98  {
99  while(!isRoot[root])
100  root++;
101  int numInAggFound = 0;
102  for(int j = 0; j < numFineNodes; j++)
103  {
104  if(j == root) //don't make a connection from the root to itself
105  {
106  numInAggFound++;
107  continue;
108  }
109  if(vertex2AggId[root] == vertex2AggId[j])
110  {
111  vertices.push_back(root);
112  vertices.push_back(j);
113  geomSizes.push_back(2);
114  //Also draw the free endpoint explicitly for the current line
115  vertices.push_back(j);
116  geomSizes.push_back(1);
117  numInAggFound++;
118  //if(numInAggFound == aggSizes_[vertex2AggId_[root]]) //don't spend more time looking if done with that root
119  // break;
120  }
121  }
122  root++; //get set up to look for the next root
123  }
124  }
125 
126  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
127  void VisualizationHelpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::doConvexHulls2D(std::vector<int>& vertices, std::vector<int>& geomSizes, LO numLocalAggs, LO numFineNodes, const std::vector<bool>& isRoot, const ArrayRCP<LO>& vertex2AggId, const Teuchos::ArrayRCP<const double>& xCoords, const Teuchos::ArrayRCP<const double>& yCoords, const Teuchos::ArrayRCP<const double>& zCoords) {
128  for(int agg = 0; agg < numLocalAggs; agg++) {
129  std::list<int> aggNodes;
130  for(int i = 0; i < numFineNodes; i++) {
131  if(vertex2AggId[i] == agg)
132  aggNodes.push_back(i);
133  }
134  //have a list of nodes in the aggregate
135  TEUCHOS_TEST_FOR_EXCEPTION(aggNodes.size() == 0, Exceptions::RuntimeError,
136  "CoarseningVisualization::doConvexHulls2D: aggregate contains zero nodes!");
137  if(aggNodes.size() == 1) {
138  vertices.push_back(aggNodes.front());
139  geomSizes.push_back(1);
140  continue;
141  }
142  if(aggNodes.size() == 2) {
143  vertices.push_back(aggNodes.front());
144  vertices.push_back(aggNodes.back());
145  geomSizes.push_back(2);
146  continue;
147  }
148  //check if all points are collinear, need to explicitly draw a line in that case.
149  bool collinear = true; //assume true at first, if a segment not parallel to others then clear
150  {
151  std::list<int>::iterator it = aggNodes.begin();
152  myVec3 firstPoint(xCoords[*it], yCoords[*it], 0);
153  it++;
154  myVec3 secondPoint(xCoords[*it], yCoords[*it], 0);
155  it++; //it now points to third node in the aggregate
156  myVec3 norm1(-(secondPoint.y - firstPoint.y), secondPoint.x - firstPoint.x, 0);
157  do {
158  myVec3 thisNorm(yCoords[*it] - firstPoint.y, firstPoint.x - xCoords[*it], 0);
159  //rotate one of the vectors by 90 degrees so that dot product is 0 if the two are parallel
160  double temp = thisNorm.x;
161  thisNorm.x = thisNorm.y;
162  thisNorm.y = temp;
163  double comp = dotProduct(norm1, thisNorm);
164  if(-1e-8 > comp || comp > 1e-8) {
165  collinear = false;
166  break;
167  }
168  it++;
169  }
170  while(it != aggNodes.end());
171  }
172  if(collinear)
173  {
174  //find the most distant two points in the plane and use as endpoints of line representing agg
175  std::list<int>::iterator min = aggNodes.begin(); //min X then min Y where x is a tie
176  std::list<int>::iterator max = aggNodes.begin(); //max X then max Y where x is a tie
177  for(std::list<int>::iterator it = ++aggNodes.begin(); it != aggNodes.end(); it++) {
178  if(xCoords[*it] < xCoords[*min])
179  min = it;
180  else if(xCoords[*it] == xCoords[*min]) {
181  if(yCoords[*it] < yCoords[*min])
182  min = it;
183  }
184  if(xCoords[*it] > xCoords[*max])
185  max = it;
186  else if(xCoords[*it] == xCoords[*max]) {
187  if(yCoords[*it] > yCoords[*max])
188  max = it;
189  }
190  }
191  //Just set up a line between nodes *min and *max
192  vertices.push_back(*min);
193  vertices.push_back(*max);
194  geomSizes.push_back(2);
195  continue; //jump to next aggregate in loop
196  }
197  std::vector<myVec2> points;
198  std::vector<int> nodes;
199  for(std::list<int>::iterator it = aggNodes.begin(); it != aggNodes.end(); it++) {
200  points.push_back(myVec2(xCoords[*it], yCoords[*it]));
201  nodes.push_back(*it);
202  }
203  std::vector<int> hull = giftWrap(points, nodes, xCoords, yCoords);
204  vertices.reserve(vertices.size() + hull.size());
205  for(size_t i = 0; i < hull.size(); i++) {
206  vertices.push_back(hull[i]);
207  }
208  geomSizes.push_back(hull.size());
209  }
210  }
211 
212  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
213  void VisualizationHelpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::doConvexHulls3D(std::vector<int>& vertices, std::vector<int>& geomSizes, LO numLocalAggs, LO numFineNodes, const std::vector<bool>& isRoot, const ArrayRCP<LO>& vertex2AggId, const Teuchos::ArrayRCP<const double>& xCoords, const Teuchos::ArrayRCP<const double>& yCoords, const Teuchos::ArrayRCP<const double>& zCoords) {
214  //Use 3D quickhull algo.
215  //Vector of node indices representing triangle vertices
216  //Note: Calculate the hulls first since will only include point data for points in the hulls
217  //Effectively the size() of vertIndices after each hull is added to it
218  typedef std::list<int>::iterator Iter;
219  for(int agg = 0; agg < numLocalAggs; agg++) {
220  std::list<int> aggNodes; //At first, list of all nodes in the aggregate. As nodes are enclosed or included by/in hull, remove them
221  for(int i = 0; i < numFineNodes; i++) {
222  if(vertex2AggId[i] == agg)
223  aggNodes.push_back(i);
224  }
225  //First, check anomalous cases
226  TEUCHOS_TEST_FOR_EXCEPTION(aggNodes.size() == 0, Exceptions::RuntimeError,
227  "CoarseningVisualization::doConvexHulls3D: aggregate contains zero nodes!");
228  if(aggNodes.size() == 1) {
229  vertices.push_back(aggNodes.front());
230  geomSizes.push_back(1);
231  continue;
232  } else if(aggNodes.size() == 2) {
233  vertices.push_back(aggNodes.front());
234  vertices.push_back(aggNodes.back());
235  geomSizes.push_back(2);
236  continue;
237  }
238  //check for collinearity
239  bool areCollinear = true;
240  {
241  Iter it = aggNodes.begin();
242  myVec3 firstVec(xCoords[*it], yCoords[*it], zCoords[*it]);
243  myVec3 comp;
244  {
245  it++;
246  myVec3 secondVec(xCoords[*it], yCoords[*it], zCoords[*it]); //cross this with other vectors to compare
247  comp = vecSubtract(secondVec, firstVec);
248  it++;
249  }
250  for(; it != aggNodes.end(); it++) {
251  myVec3 thisVec(xCoords[*it], yCoords[*it], zCoords[*it]);
252  myVec3 cross = crossProduct(vecSubtract(thisVec, firstVec), comp);
253  if(mymagnitude(cross) > 1e-10) {
254  areCollinear = false;
255  break;
256  }
257  }
258  }
259  if(areCollinear)
260  {
261  //find the endpoints of segment describing all the points
262  //compare x, if tie compare y, if tie compare z
263  Iter min = aggNodes.begin();
264  Iter max = aggNodes.begin();
265  Iter it = ++aggNodes.begin();
266  for(; it != aggNodes.end(); it++) {
267  if(xCoords[*it] < xCoords[*min]) min = it;
268  else if(xCoords[*it] == xCoords[*min]) {
269  if(yCoords[*it] < yCoords[*min]) min = it;
270  else if(yCoords[*it] == yCoords[*min]) {
271  if(zCoords[*it] < zCoords[*min]) min = it;
272  }
273  }
274  if(xCoords[*it] > xCoords[*max]) max = it;
275  else if(xCoords[*it] == xCoords[*max]) {
276  if(yCoords[*it] > yCoords[*max]) max = it;
277  else if(yCoords[*it] == yCoords[*max]) {
278  if(zCoords[*it] > zCoords[*max])
279  max = it;
280  }
281  }
282  }
283  vertices.push_back(*min);
284  vertices.push_back(*max);
285  geomSizes.push_back(2);
286  continue;
287  }
288  bool areCoplanar = true;
289  {
290  //number of points is known to be >= 3
291  Iter vert = aggNodes.begin();
292  myVec3 v1(xCoords[*vert], yCoords[*vert], zCoords[*vert]);
293  vert++;
294  myVec3 v2(xCoords[*vert], yCoords[*vert], zCoords[*vert]);
295  vert++;
296  myVec3 v3(xCoords[*vert], yCoords[*vert], zCoords[*vert]);
297  vert++;
298  //Make sure the first three points aren't also collinear (need a non-degenerate triangle to get a normal)
299  while(mymagnitude(crossProduct(vecSubtract(v1, v2), vecSubtract(v1, v3))) < 1e-10) {
300  //Replace the third point with the next point
301  v3 = myVec3(xCoords[*vert], yCoords[*vert], zCoords[*vert]);
302  vert++;
303  }
304  for(; vert != aggNodes.end(); vert++) {
305  myVec3 pt(xCoords[*vert], yCoords[*vert], zCoords[*vert]);
306  if(fabs(pointDistFromTri(pt, v1, v2, v3)) > 1e-12) {
307  areCoplanar = false;
308  break;
309  }
310  }
311  if(areCoplanar) {
312  //do 2D convex hull
313  myVec3 planeNorm = getNorm(v1, v2, v3);
314  planeNorm.x = fabs(planeNorm.x);
315  planeNorm.y = fabs(planeNorm.y);
316  planeNorm.z = fabs(planeNorm.z);
317  std::vector<myVec2> points;
318  std::vector<int> nodes;
319  if(planeNorm.x >= planeNorm.y && planeNorm.x >= planeNorm.z) {
320  //project points to yz plane to make hull
321  for(Iter it = aggNodes.begin(); it != aggNodes.end(); it++) {
322  nodes.push_back(*it);
323  points.push_back(myVec2(yCoords[*it], zCoords[*it]));
324  }
325  }
326  if(planeNorm.y >= planeNorm.x && planeNorm.y >= planeNorm.z) {
327  //use xz
328  for(Iter it = aggNodes.begin(); it != aggNodes.end(); it++) {
329  nodes.push_back(*it);
330  points.push_back(myVec2(xCoords[*it], zCoords[*it]));
331  }
332  }
333  if(planeNorm.z >= planeNorm.x && planeNorm.z >= planeNorm.y) {
334  for(Iter it = aggNodes.begin(); it != aggNodes.end(); it++) {
335  nodes.push_back(*it);
336  points.push_back(myVec2(xCoords[*it], yCoords[*it]));
337  }
338  }
339  std::vector<int> convhull2d = giftWrap(points, nodes, xCoords, yCoords);
340  geomSizes.push_back(convhull2d.size());
341  vertices.reserve(vertices.size() + convhull2d.size());
342  for(size_t i = 0; i < convhull2d.size(); i++)
343  vertices.push_back(convhull2d[i]);
344  continue;
345  }
346  }
347  Iter exIt = aggNodes.begin(); //iterator to be used for searching for min/max x/y/z
348  int extremeSix[] = {*exIt, *exIt, *exIt, *exIt, *exIt, *exIt}; //nodes with minimumX, maxX, minY ...
349  exIt++;
350  for(; exIt != aggNodes.end(); exIt++) {
351  Iter it = exIt;
352  if(xCoords[*it] < xCoords[extremeSix[0]] ||
353  (xCoords[*it] == xCoords[extremeSix[0]] && yCoords[*it] < yCoords[extremeSix[0]]) ||
354  (xCoords[*it] == xCoords[extremeSix[0]] && yCoords[*it] == yCoords[extremeSix[0]] && zCoords[*it] < zCoords[extremeSix[0]]))
355  extremeSix[0] = *it;
356  if(xCoords[*it] > xCoords[extremeSix[1]] ||
357  (xCoords[*it] == xCoords[extremeSix[1]] && yCoords[*it] > yCoords[extremeSix[1]]) ||
358  (xCoords[*it] == xCoords[extremeSix[1]] && yCoords[*it] == yCoords[extremeSix[1]] && zCoords[*it] > zCoords[extremeSix[1]]))
359  extremeSix[1] = *it;
360  if(yCoords[*it] < yCoords[extremeSix[2]] ||
361  (yCoords[*it] == yCoords[extremeSix[2]] && zCoords[*it] < zCoords[extremeSix[2]]) ||
362  (yCoords[*it] == yCoords[extremeSix[2]] && zCoords[*it] == zCoords[extremeSix[2]] && xCoords[*it] < xCoords[extremeSix[2]]))
363  extremeSix[2] = *it;
364  if(yCoords[*it] > yCoords[extremeSix[3]] ||
365  (yCoords[*it] == yCoords[extremeSix[3]] && zCoords[*it] > zCoords[extremeSix[3]]) ||
366  (yCoords[*it] == yCoords[extremeSix[3]] && zCoords[*it] == zCoords[extremeSix[3]] && xCoords[*it] > xCoords[extremeSix[3]]))
367  extremeSix[3] = *it;
368  if(zCoords[*it] < zCoords[extremeSix[4]] ||
369  (zCoords[*it] == zCoords[extremeSix[4]] && xCoords[*it] < xCoords[extremeSix[4]]) ||
370  (zCoords[*it] == zCoords[extremeSix[4]] && xCoords[*it] == xCoords[extremeSix[4]] && yCoords[*it] < yCoords[extremeSix[4]]))
371  extremeSix[4] = *it;
372  if(zCoords[*it] > zCoords[extremeSix[5]] ||
373  (zCoords[*it] == zCoords[extremeSix[5]] && xCoords[*it] > xCoords[extremeSix[5]]) ||
374  (zCoords[*it] == zCoords[extremeSix[5]] && xCoords[*it] == xCoords[extremeSix[5]] && yCoords[*it] > zCoords[extremeSix[5]]))
375  extremeSix[5] = *it;
376  }
377  myVec3 extremeVectors[6];
378  for(int i = 0; i < 6; i++) {
379  myVec3 thisExtremeVec(xCoords[extremeSix[i]], yCoords[extremeSix[i]], zCoords[extremeSix[i]]);
380  extremeVectors[i] = thisExtremeVec;
381  }
382  double maxDist = 0;
383  int base1 = 0; //ints from 0-5: which pair out of the 6 extreme points are the most distant? (indices in extremeSix and extremeVectors)
384  int base2 = 0;
385  for(int i = 0; i < 5; i++) {
386  for(int j = i + 1; j < 6; j++) {
387  double thisDist = distance(extremeVectors[i], extremeVectors[j]);
388  if(thisDist > maxDist) {
389  maxDist = thisDist;
390  base1 = i;
391  base2 = j;
392  }
393  }
394  }
395  std::list<myTriangle> hullBuilding; //each Triangle is a triplet of nodes (int IDs) that form a triangle
396  //remove base1 and base2 iters from aggNodes, they are known to be in the hull
397  aggNodes.remove(extremeSix[base1]);
398  aggNodes.remove(extremeSix[base2]);
399  //extremeSix[base1] and [base2] still have the myVec3 representation
400  myTriangle tri;
401  tri.v1 = extremeSix[base1];
402  tri.v2 = extremeSix[base2];
403  //Now find the point that is furthest away from the line between base1 and base2
404  maxDist = 0;
405  //need the vectors to do "quadruple product" formula
406  myVec3 b1 = extremeVectors[base1];
407  myVec3 b2 = extremeVectors[base2];
408  Iter thirdNode;
409  for(Iter node = aggNodes.begin(); node != aggNodes.end(); node++) {
410  myVec3 nodePos(xCoords[*node], yCoords[*node], zCoords[*node]);
411  double dist = mymagnitude(crossProduct(vecSubtract(nodePos, b1), vecSubtract(nodePos, b2))) / mymagnitude(vecSubtract(b2, b1));
412  if(dist > maxDist) {
413  maxDist = dist;
414  thirdNode = node;
415  }
416  }
417  //Now know the last node in the first triangle
418  tri.v3 = *thirdNode;
419  hullBuilding.push_back(tri);
420  myVec3 b3(xCoords[*thirdNode], yCoords[*thirdNode], zCoords[*thirdNode]);
421  aggNodes.erase(thirdNode);
422  //Find the fourth node (most distant from triangle) to form tetrahedron
423  maxDist = 0;
424  int fourthVertex = -1;
425  for(Iter node = aggNodes.begin(); node != aggNodes.end(); node++) {
426  myVec3 thisNode(xCoords[*node], yCoords[*node], zCoords[*node]);
427  double nodeDist = pointDistFromTri(thisNode, b1, b2, b3);
428  if(nodeDist > maxDist) {
429  maxDist = nodeDist;
430  fourthVertex = *node;
431  }
432  }
433  aggNodes.remove(fourthVertex);
434  myVec3 b4(xCoords[fourthVertex], yCoords[fourthVertex], zCoords[fourthVertex]);
435  //Add three new triangles to hullBuilding to form the first tetrahedron
436  //use tri to hold the triangle info temporarily before being added to list
437  tri = hullBuilding.front();
438  tri.v1 = fourthVertex;
439  hullBuilding.push_back(tri);
440  tri = hullBuilding.front();
441  tri.v2 = fourthVertex;
442  hullBuilding.push_back(tri);
443  tri = hullBuilding.front();
444  tri.v3 = fourthVertex;
445  hullBuilding.push_back(tri);
446  //now orient all four triangles so that the vertices are oriented clockwise (so getNorm_ points outward for each)
447  myVec3 barycenter((b1.x + b2.x + b3.x + b4.x) / 4.0, (b1.y + b2.y + b3.y + b4.y) / 4.0, (b1.z + b2.z + b3.z + b4.z) / 4.0);
448  for(std::list<myTriangle>::iterator tetTri = hullBuilding.begin(); tetTri != hullBuilding.end(); tetTri++) {
449  myVec3 triVert1(xCoords[tetTri->v1], yCoords[tetTri->v1], zCoords[tetTri->v1]);
450  myVec3 triVert2(xCoords[tetTri->v2], yCoords[tetTri->v2], zCoords[tetTri->v2]);
451  myVec3 triVert3(xCoords[tetTri->v3], yCoords[tetTri->v3], zCoords[tetTri->v3]);
452  myVec3 trinorm = getNorm(triVert1, triVert2, triVert3);
453  myVec3 ptInPlane = tetTri == hullBuilding.begin() ? b1 : b4; //first triangle definitely has b1 in it, other three definitely have b4
454  if(isInFront(barycenter, ptInPlane, trinorm)) {
455  //don't want the faces of the tetrahedron to face inwards (towards barycenter) so reverse orientation
456  //by swapping two vertices
457  int temp = tetTri->v1;
458  tetTri->v1 = tetTri->v2;
459  tetTri->v2 = temp;
460  }
461  }
462  //now, have starting polyhedron in hullBuilding (all faces are facing outwards according to getNorm_) and remaining nodes to process are in aggNodes
463  //recursively, for each triangle, make a list of the points that are 'in front' of the triangle. Find the point with the maximum distance from the triangle.
464  //Add three new triangles, each sharing one edge with the original triangle but now with the most distant point as a vertex. Remove the most distant point from
465  //the list of all points that need to be processed. Also from that list remove all points that are in front of the original triangle but not in front of all three
466  //new triangles, since they are now enclosed in the hull.
467  //Construct point lists for each face of the tetrahedron individually.
468  myVec3 trinorms[4]; //normals to the four tetrahedron faces, now oriented outwards
469  int index = 0;
470  for(std::list<myTriangle>::iterator tetTri = hullBuilding.begin(); tetTri != hullBuilding.end(); tetTri++) {
471  myVec3 triVert1(xCoords[tetTri->v1], yCoords[tetTri->v1], zCoords[tetTri->v1]);
472  myVec3 triVert2(xCoords[tetTri->v2], yCoords[tetTri->v2], zCoords[tetTri->v2]);
473  myVec3 triVert3(xCoords[tetTri->v3], yCoords[tetTri->v3], zCoords[tetTri->v3]);
474  trinorms[index] = getNorm(triVert1, triVert2, triVert3);
475  index++;
476  }
477  std::list<int> startPoints1;
478  std::list<int> startPoints2;
479  std::list<int> startPoints3;
480  std::list<int> startPoints4;
481  //scope this so that 'point' is not in function scope
482  {
483  Iter point = aggNodes.begin();
484  while(!aggNodes.empty()) //this removes points one at a time as they are put in startPointsN or are already done
485  {
486  myVec3 pointVec(xCoords[*point], yCoords[*point], zCoords[*point]);
487  //Note: Because of the way the tetrahedron faces are constructed above,
488  //face 1 definitely contains b1 and faces 2-4 definitely contain b4.
489  if(isInFront(pointVec, b1, trinorms[0])) {
490  startPoints1.push_back(*point);
491  point = aggNodes.erase(point);
492  } else if(isInFront(pointVec, b4, trinorms[1])) {
493  startPoints2.push_back(*point);
494  point = aggNodes.erase(point);
495  } else if(isInFront(pointVec, b4, trinorms[2])) {
496  startPoints3.push_back(*point);
497  point = aggNodes.erase(point);
498  } else if(isInFront(pointVec, b4, trinorms[3])) {
499  startPoints4.push_back(*point);
500  point = aggNodes.erase(point);
501  } else {
502  point = aggNodes.erase(point); //points here are already inside tetrahedron.
503  }
504  }
505  //Call processTriangle for each triangle in the initial tetrahedron, one at a time.
506  }
507  typedef std::list<myTriangle>::iterator TriIter;
508  TriIter firstTri = hullBuilding.begin();
509  myTriangle start1 = *firstTri;
510  firstTri++;
511  myTriangle start2 = *firstTri;
512  firstTri++;
513  myTriangle start3 = *firstTri;
514  firstTri++;
515  myTriangle start4 = *firstTri;
516  //kick off depth-first recursive filling of hullBuilding list with all triangles in the convex hull
517  if(!startPoints1.empty())
518  processTriangle(hullBuilding, start1, startPoints1, barycenter, xCoords, yCoords, zCoords);
519  if(!startPoints2.empty())
520  processTriangle(hullBuilding, start2, startPoints2, barycenter, xCoords, yCoords, zCoords);
521  if(!startPoints3.empty())
522  processTriangle(hullBuilding, start3, startPoints3, barycenter, xCoords, yCoords, zCoords);
523  if(!startPoints4.empty())
524  processTriangle(hullBuilding, start4, startPoints4, barycenter, xCoords, yCoords, zCoords);
525  //hullBuilding now has all triangles that make up this hull.
526  //Dump hullBuilding info into the list of all triangles for the scene.
527  vertices.reserve(vertices.size() + 3 * hullBuilding.size());
528  for(TriIter hullTri = hullBuilding.begin(); hullTri != hullBuilding.end(); hullTri++) {
529  vertices.push_back(hullTri->v1);
530  vertices.push_back(hullTri->v2);
531  vertices.push_back(hullTri->v3);
532  geomSizes.push_back(3);
533  }
534  }
535  }
536 
537  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
538  void VisualizationHelpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::doGraphEdges(std::vector<int>& vertices, std::vector<int>& geomSizes, Teuchos::RCP<GraphBase>& G, Teuchos::ArrayRCP<const double> & fx, Teuchos::ArrayRCP<const double> & fy, Teuchos::ArrayRCP<const double> & fz) {
539  ArrayView<const Scalar> values;
540  ArrayView<const LocalOrdinal> neighbors;
541 
542  std::vector<std::pair<int, int> > vert1; //vertices (node indices)
543 
544  ArrayView<const LocalOrdinal> indices;
545  for(LocalOrdinal locRow = 0; locRow < LocalOrdinal(G->GetNodeNumVertices()); locRow++) {
546  neighbors = G->getNeighborVertices(locRow);
547  //Add those local indices (columns) to the list of connections (which are pairs of the form (localM, localN))
548  for(int gEdge = 0; gEdge < int(neighbors.size()); ++gEdge) {
549  vert1.push_back(std::pair<int, int>(locRow, neighbors[gEdge]));
550  }
551  }
552  for(size_t i = 0; i < vert1.size(); i ++) {
553  if(vert1[i].first > vert1[i].second) {
554  int temp = vert1[i].first;
555  vert1[i].first = vert1[i].second;
556  vert1[i].second = temp;
557  }
558  }
559  std::sort(vert1.begin(), vert1.end());
560  std::vector<std::pair<int, int> >::iterator newEnd = unique(vert1.begin(), vert1.end()); //remove duplicate edges
561  vert1.erase(newEnd, vert1.end());
562  //std::vector<int> points1;
563  vertices.reserve(2 * vert1.size());
564  geomSizes.reserve(vert1.size());
565  for(size_t i = 0; i < vert1.size(); i++) {
566  vertices.push_back(vert1[i].first);
567  vertices.push_back(vert1[i].second);
568  geomSizes.push_back(2);
569  }
570  }
571 
572  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
574  {
575  return myVec3(v1.y * v2.z - v1.z * v2.y, v1.z * v2.x - v1.x * v2.z, v1.x * v2.y - v1.y * v2.x);
576  }
577 
578  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
580  {
581  return v1.x * v2.x + v1.y * v2.y;
582  }
583 
584  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
586  {
587  return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z;
588  }
589 
590  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
592  {
593  myVec3 rel(point.x - inPlane.x, point.y - inPlane.y, point.z - inPlane.z); //position of the point relative to the plane
594  return dotProduct(rel, n) > 1e-12 ? true : false;
595  }
596 
597  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
599  {
600  return sqrt(dotProduct(vec, vec));
601  }
602 
603  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
605  {
606  return sqrt(dotProduct(vec, vec));
607  }
608 
609  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
611  {
612  return mymagnitude(vecSubtract(p1, p2));
613  }
614 
615 
616  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
618  {
619  return mymagnitude(vecSubtract(p1, p2));
620  }
621 
622  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
624  {
625  return myVec2(v1.x - v2.x, v1.y - v2.y);
626  }
627 
628  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
630  {
631  return myVec3(v1.x - v2.x, v1.y - v2.y, v1.z - v2.z);
632  }
633 
634  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
635  myVec2 VisualizationHelpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getNorm(myVec2 v) //"normal" to a 2D vector - just rotate 90 degrees to left
636  {
637  return myVec2(v.y, -v.x);
638  }
639 
640  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
641  myVec3 VisualizationHelpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getNorm(myVec3 v1, myVec3 v2, myVec3 v3) //normal to face of triangle (will be outward rel. to polyhedron) (v1, v2, v3 are in CCW order when normal is toward viewpoint)
642  {
643  return crossProduct(vecSubtract(v2, v1), vecSubtract(v3, v1));
644  }
645 
646  //get minimum distance from 'point' to plane containing v1, v2, v3 (or the triangle with v1, v2, v3 as vertices)
647  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
649  {
650  using namespace std;
651  myVec3 norm = getNorm(v1, v2, v3);
652  //must normalize the normal vector
653  double normScl = mymagnitude(norm);
654  double rv = 0.0;
655  if (normScl > 1e-8) {
656  norm.x /= normScl;
657  norm.y /= normScl;
658  norm.z /= normScl;
659  rv = fabs(dotProduct(norm, vecSubtract(point, v1)));
660  } else {
661  // triangle is degenerated
662  myVec3 test1 = vecSubtract(v3, v1);
663  myVec3 test2 = vecSubtract(v2, v1);
664  bool useTest1 = mymagnitude(test1) > 0.0 ? true : false;
665  bool useTest2 = mymagnitude(test2) > 0.0 ? true : false;
666  if(useTest1 == true) {
667  double normScl1 = mymagnitude(test1);
668  test1.x /= normScl1;
669  test1.y /= normScl1;
670  test1.z /= normScl1;
671  rv = fabs(dotProduct(test1, vecSubtract(point, v1)));
672  } else if (useTest2 == true) {
673  double normScl2 = mymagnitude(test2);
674  test2.x /= normScl2;
675  test2.y /= normScl2;
676  test2.z /= normScl2;
677  rv = fabs(dotProduct(test2, vecSubtract(point, v1)));
678  } else {
679  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError,
680  "VisualizationHelpers::pointDistFromTri: Could not determine the distance of a point to a triangle.");
681  }
682 
683  }
684  return rv;
685  }
686 
687  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
688  std::vector<myTriangle> VisualizationHelpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::processTriangle(std::list<myTriangle>& tris, myTriangle tri, std::list<int>& pointsInFront, myVec3& barycenter, const Teuchos::ArrayRCP<const double>& xCoords, const Teuchos::ArrayRCP<const double>& yCoords, const Teuchos::ArrayRCP<const double>& zCoords) {
689  //*tri is in the tris list, and is the triangle to process here. tris is a complete list of all triangles in the hull so far. pointsInFront is only a list of the nodes in front of tri. Need coords also.
690  //precondition: each triangle is already oriented so that getNorm_(v1, v2, v3) points outward (away from interior of hull)
691  //First find the point furthest from triangle.
692  using namespace std;
693  typedef std::list<int>::iterator Iter;
694  typedef std::list<myTriangle>::iterator TriIter;
695  typedef list<pair<int, int> >::iterator EdgeIter;
696  double maxDist = 0;
697  //Need vector representations of triangle's vertices
698  myVec3 v1(xCoords[tri.v1], yCoords[tri.v1], zCoords[tri.v1]);
699  myVec3 v2(xCoords[tri.v2], yCoords[tri.v2], zCoords[tri.v2]);
700  myVec3 v3(xCoords[tri.v3], yCoords[tri.v3], zCoords[tri.v3]);
701  myVec3 farPointVec; //useful to have both the point's coordinates and it's position in the list
702  Iter farPoint = pointsInFront.begin();
703  for(Iter point = pointsInFront.begin(); point != pointsInFront.end(); point++)
704  {
705  myVec3 pointVec(xCoords[*point], yCoords[*point], zCoords[*point]);
706  double dist = pointDistFromTri(pointVec, v1, v2, v3);
707  if(dist > maxDist)
708  {
709  dist = maxDist;
710  farPointVec = pointVec;
711  farPoint = point;
712  }
713  }
714  //Find all the triangles that the point is in front of (can be more than 1)
715  //At the same time, remove them from tris, as every one will be replaced later
716  vector<myTriangle> visible; //use a list of iterators so that the underlying object is still in tris
717  for(TriIter it = tris.begin(); it != tris.end();)
718  {
719  myVec3 vec1(xCoords[it->v1], yCoords[it->v1], zCoords[it->v1]);
720  myVec3 vec2(xCoords[it->v2], yCoords[it->v2], zCoords[it->v2]);
721  myVec3 vec3(xCoords[it->v3], yCoords[it->v3], zCoords[it->v3]);
722  myVec3 norm = getNorm(vec1, vec2, vec3);
723  if(isInFront(farPointVec, vec1, norm))
724  {
725  visible.push_back(*it);
726  it = tris.erase(it);
727  }
728  else
729  it++;
730  }
731  //Figure out what triangles need to be destroyed/created
732  //First create a list of edges (as std::pair<int, int>, where the two ints are the node endpoints)
733  list<pair<int, int> > horizon;
734  //For each triangle, add edges to the list iff the edge only appears once in the set of all
735  //Have members of horizon have the lower node # first, and the higher one second
736  for(vector<myTriangle>::iterator it = visible.begin(); it != visible.end(); it++)
737  {
738  pair<int, int> e1(it->v1, it->v2);
739  pair<int, int> e2(it->v2, it->v3);
740  pair<int, int> e3(it->v1, it->v3);
741  //"sort" the pair values
742  if(e1.first > e1.second)
743  {
744  int temp = e1.first;
745  e1.first = e1.second;
746  e1.second = temp;
747  }
748  if(e2.first > e2.second)
749  {
750  int temp = e2.first;
751  e2.first = e2.second;
752  e2.second = temp;
753  }
754  if(e3.first > e3.second)
755  {
756  int temp = e3.first;
757  e3.first = e3.second;
758  e3.second = temp;
759  }
760  horizon.push_back(e1);
761  horizon.push_back(e2);
762  horizon.push_back(e3);
763  }
764  //sort based on lower node first, then higher node (lexicographical ordering built in to pair)
765  horizon.sort();
766  //Remove all edges from horizon, except those that appear exactly once
767  {
768  EdgeIter it = horizon.begin();
769  while(it != horizon.end())
770  {
771  int occur = count(horizon.begin(), horizon.end(), *it);
772  if(occur > 1)
773  {
774  pair<int, int> removeVal = *it;
775  while(removeVal == *it && !(it == horizon.end()))
776  it = horizon.erase(it);
777  }
778  else
779  it++;
780  }
781  }
782  //Now make a list of new triangles being created, each of which take 2 vertices from an edge and one from farPoint
783  list<myTriangle> newTris;
784  for(EdgeIter it = horizon.begin(); it != horizon.end(); it++)
785  {
786  myTriangle t(it->first, it->second, *farPoint);
787  newTris.push_back(t);
788  }
789  //Ensure every new triangle is oriented outwards, using the barycenter of the initial tetrahedron
790  vector<myTriangle> trisToProcess;
791  vector<list<int> > newFrontPoints;
792  for(TriIter it = newTris.begin(); it != newTris.end(); it++)
793  {
794  myVec3 t1(xCoords[it->v1], yCoords[it->v1], zCoords[it->v1]);
795  myVec3 t2(xCoords[it->v2], yCoords[it->v2], zCoords[it->v2]);
796  myVec3 t3(xCoords[it->v3], yCoords[it->v3], zCoords[it->v3]);
797  if(isInFront(barycenter, t1, getNorm(t1, t2, t3)))
798  {
799  //need to swap two vertices to flip orientation of triangle
800  int temp = it->v1;
801  myVec3 tempVec = t1;
802  it->v1 = it->v2;
803  t1 = t2;
804  it->v2 = temp;
805  t2 = tempVec;
806  }
807  myVec3 outwardNorm = getNorm(t1, t2, t3); //now definitely points outwards
808  //Add the triangle to tris
809  tris.push_back(*it);
810  trisToProcess.push_back(tris.back());
811  //Make a list of the points that are in front of nextToProcess, to be passed in for processing
812  list<int> newInFront;
813  for(Iter point = pointsInFront.begin(); point != pointsInFront.end();)
814  {
815  myVec3 pointVec(xCoords[*point], yCoords[*point], zCoords[*point]);
816  if(isInFront(pointVec, t1, outwardNorm))
817  {
818  newInFront.push_back(*point);
819  point = pointsInFront.erase(point);
820  }
821  else
822  point++;
823  }
824  newFrontPoints.push_back(newInFront);
825  }
826  vector<myTriangle> allRemoved; //list of all invalid iterators that were erased by calls to processmyTriangle below
827  for(int i = 0; i < int(trisToProcess.size()); i++)
828  {
829  if(!newFrontPoints[i].empty())
830  {
831  //Comparing the 'triangle to process' to the one for this call prevents infinite recursion/stack overflow.
832  //TODO: Why was it doing that? Rounding error? Make more robust fix. But this does work for the time being.
833  if(find(allRemoved.begin(), allRemoved.end(), trisToProcess[i]) == allRemoved.end() && !(trisToProcess[i] == tri))
834  {
835  vector<myTriangle> removedList = processTriangle(tris, trisToProcess[i], newFrontPoints[i], barycenter, xCoords, yCoords, zCoords);
836  for(int j = 0; j < int(removedList.size()); j++)
837  allRemoved.push_back(removedList[j]);
838  }
839  }
840  }
841  return visible;
842  }
843 
844  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
845  std::vector<int> VisualizationHelpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::giftWrap(std::vector<myVec2>& points, std::vector<int>& nodes, const Teuchos::ArrayRCP<const double> & xCoords, const Teuchos::ArrayRCP<const double> & yCoords) {
846  TEUCHOS_TEST_FOR_EXCEPTION(points.size() < 3, Exceptions::RuntimeError,
847  "CoarseningVisualization::giftWrap: Gift wrap algorithm input has to have at least 3 points!");
848 
849 #if 1 // TAW's version to determine "minimal" node
850  // determine minimal x and y coordinates
851  double min_x =points[0].x;
852  double min_y =points[0].y;
853  for(std::vector<int>::iterator it = nodes.begin(); it != nodes.end(); it++) {
854  int i = it - nodes.begin();
855  if(points[i].x < min_x) min_x = points[i].x;
856  if(points[i].y < min_y) min_y = points[i].y;
857  }
858  // create dummy min coordinates
859  min_x -= 1.0;
860  min_y -= 1.0;
861  myVec2 dummy_min(min_x, min_y);
862 
863  // loop over all nodes and determine nodes with minimal distance to (min_x, min_y)
864  std::vector<int> hull;
865  myVec2 min = points[0];
866  double mindist = distance(min,dummy_min);
867  std::vector<int>::iterator minNode = nodes.begin();
868  for(std::vector<int>::iterator it = nodes.begin(); it != nodes.end(); it++) {
869  int i = it - nodes.begin();
870  if(distance(points[i],dummy_min) < mindist) {
871  mindist = distance(points[i],dummy_min);
872  min = points[i];
873  minNode = it;
874  }
875  }
876  hull.push_back(*minNode);
877 #else // Brian's code
878  std::vector<int> hull;
879  std::vector<int>::iterator minNode = nodes.begin();
880  myVec2 min = points[0];
881  for(std::vector<int>::iterator it = nodes.begin(); it != nodes.end(); it++)
882  {
883  int i = it - nodes.begin();
884  if(points[i].x < min.x || (fabs(points[i].x - min.x) < 1e-10 && points[i].y < min.y))
885  {
886  min = points[i];
887  minNode = it;
888  }
889  }
890  hull.push_back(*minNode);
891 #endif
892 
893  bool includeMin = false;
894  //int debug_it = 0;
895  while(1)
896  {
897  std::vector<int>::iterator leftMost = nodes.begin();
898  if(!includeMin && leftMost == minNode)
899  {
900  leftMost++;
901  }
902  std::vector<int>::iterator it = leftMost;
903  it++;
904  for(; it != nodes.end(); it++)
905  {
906  if(it == minNode && !includeMin) //don't compare to min on very first sweep
907  continue;
908  if(*it == hull.back())
909  continue;
910  //see if it is in front of line containing nodes thisHull.back() and leftMost
911  //first get the left normal of leftMost - thisHull.back() (<dy, -dx>)
912  myVec2 leftMostVec = points[leftMost - nodes.begin()];
913  myVec2 lastVec(xCoords[hull.back()], yCoords[hull.back()]);
914  myVec2 testNorm = getNorm(vecSubtract(leftMostVec, lastVec));
915  //now dot testNorm with *it - leftMost. If dot is positive, leftMost becomes it. If dot is zero, take one further from thisHull.back().
916  myVec2 itVec(xCoords[*it], yCoords[*it]);
917  double dotProd = dotProduct(testNorm, vecSubtract(itVec, lastVec));
918  if(-1e-8 < dotProd && dotProd < 1e-8)
919  {
920  //thisHull.back(), it and leftMost are collinear.
921  //Just sum the differences in x and differences in y for each and compare to get further one, don't need distance formula
922  myVec2 itDist = vecSubtract(itVec, lastVec);
923  myVec2 leftMostDist = vecSubtract(leftMostVec, lastVec);
924  if(fabs(itDist.x) + fabs(itDist.y) > fabs(leftMostDist.x) + fabs(leftMostDist.y)) {
925  leftMost = it;
926  }
927  }
928  else if(dotProd > 0) {
929  leftMost = it;
930 
931  }
932  }
933  //if leftMost is min, then the loop is complete.
934  if(*leftMost == *minNode)
935  break;
936  hull.push_back(*leftMost);
937  includeMin = true; //have found second point (the one after min) so now include min in the searches
938  //debug_it ++;
939  //if(debug_it > 100) exit(0); //break;
940  }
941  return hull;
942  }
943 
944  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
945  std::vector<int> VisualizationHelpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::makeUnique(std::vector<int>& vertices) const
946  {
947  using namespace std;
948  vector<int> uniqueNodes = vertices;
949  sort(uniqueNodes.begin(), uniqueNodes.end());
950  vector<int>::iterator newUniqueFineEnd = unique(uniqueNodes.begin(), uniqueNodes.end());
951  uniqueNodes.erase(newUniqueFineEnd, uniqueNodes.end());
952  //uniqueNodes is now a sorted list of the nodes whose info actually goes in file
953  //Now replace values in vertices with locations of the old values in uniqueFine
954  for(int i = 0; i < int(vertices.size()); i++)
955  {
956  int lo = 0;
957  int hi = uniqueNodes.size() - 1;
958  int mid = 0;
959  int search = vertices[i];
960  while(lo <= hi)
961  {
962  mid = lo + (hi - lo) / 2;
963  if(uniqueNodes[mid] == search)
964  break;
965  else if(uniqueNodes[mid] > search)
966  hi = mid - 1;
967  else
968  lo = mid + 1;
969  }
970  if(uniqueNodes[mid] != search)
971  throw runtime_error("Issue in makeUnique_() - a point wasn't found in list.");
972  vertices[i] = mid;
973  }
974  return uniqueNodes;
975  }
976 
977  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
978  std::string VisualizationHelpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::replaceAll(std::string result, const std::string& replaceWhat, const std::string& replaceWithWhat) const {
979  while(1) {
980  const int pos = result.find(replaceWhat);
981  if (pos == -1)
982  break;
983  result.replace(pos, replaceWhat.size(), replaceWithWhat);
984  }
985  return result;
986  }
987 
988  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
989  std::string VisualizationHelpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getFileName(int numProcs, int myRank, int level, const Teuchos::ParameterList & pL) const {
990  std::string filenameToWrite = getBaseFileName(numProcs, level, pL);
991  filenameToWrite = this->replaceAll(filenameToWrite, "%PROCID", toString(myRank));
992  return filenameToWrite;
993  }
994 
995  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
996  std::string VisualizationHelpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getBaseFileName(int numProcs, int level, const Teuchos::ParameterList & pL) const {
997  std::string filenameToWrite = pL.get<std::string>("visualization: output filename");
998  int timeStep = pL.get<int>("visualization: output file: time step");
999  int iter = pL.get<int>("visualization: output file: iter");
1000 
1001  if(filenameToWrite.rfind(".vtu") == std::string::npos)
1002  filenameToWrite.append(".vtu");
1003  if(numProcs > 1 && filenameToWrite.rfind("%PROCID") == std::string::npos) //filename can't be identical between processsors in parallel problem
1004  filenameToWrite.insert(filenameToWrite.rfind(".vtu"), "-proc%PROCID");
1005 
1006  filenameToWrite = this->replaceAll(filenameToWrite, "%LEVELID", toString(level));
1007  filenameToWrite = this->replaceAll(filenameToWrite, "%TIMESTEP", toString(timeStep));
1008  filenameToWrite = this->replaceAll(filenameToWrite, "%ITER", toString(iter));
1009  return filenameToWrite;
1010  }
1011 
1012  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1013  std::string VisualizationHelpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getPVTUFileName(int numProcs, int myRank, int level, const Teuchos::ParameterList & pL) const {
1014  std::string filenameToWrite = getBaseFileName(numProcs, level, pL);
1015  std::string masterStem = filenameToWrite.substr(0, filenameToWrite.rfind(".vtu"));
1016  masterStem = this->replaceAll(masterStem, "%PROCID", "");
1017  std::string pvtuFilename = masterStem + "-master.pvtu";
1018  return pvtuFilename;
1019  }
1020 
1021  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1022  void VisualizationHelpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeFileVTKOpening(std::ofstream & fout, std::vector<int> & uniqueFine, std::vector<int> & geomSizesFine) const {
1023  std::string styleName = "PointCloud"; // TODO adapt this
1024 
1025  std::string indent = " ";
1026  fout << "<!--" << styleName << " Aggregates Visualization-->" << std::endl;
1027  fout << "<VTKFile type=\"UnstructuredGrid\" byte_order=\"LittleEndian\">" << std::endl;
1028  fout << " <UnstructuredGrid>" << std::endl;
1029  fout << " <Piece NumberOfPoints=\"" << uniqueFine.size() << "\" NumberOfCells=\"" << geomSizesFine.size() << "\">" << std::endl;
1030  fout << " <PointData Scalars=\"Node Aggregate Processor\">" << std::endl;
1031  }
1032 
1033  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1034  void VisualizationHelpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeFileVTKNodes(std::ofstream & fout, std::vector<int> & uniqueFine, Teuchos::RCP<const Map> & nodeMap) const {
1035  std::string indent = " ";
1036  fout << " <DataArray type=\"Int32\" Name=\"Node\" format=\"ascii\">" << std::endl;
1037  indent = " ";
1038  bool localIsGlobal = GlobalOrdinal(nodeMap->getGlobalNumElements()) == GlobalOrdinal(nodeMap->getNodeNumElements());
1039  for(size_t i = 0; i < uniqueFine.size(); i++)
1040  {
1041  if(localIsGlobal)
1042  {
1043  fout << uniqueFine[i] << " "; //if all nodes are on this processor, do not map from local to global
1044  }
1045  else
1046  fout << nodeMap->getGlobalElement(uniqueFine[i]) << " ";
1047  if(i % 10 == 9)
1048  fout << std::endl << indent;
1049  }
1050  fout << std::endl;
1051  fout << " </DataArray>" << std::endl;
1052  }
1053 
1054  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1055  void VisualizationHelpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeFileVTKData(std::ofstream & fout, std::vector<int> & uniqueFine, LocalOrdinal myAggOffset, ArrayRCP<LocalOrdinal> & vertex2AggId, int myRank) const {
1056  std::string indent = " ";
1057  fout << " <DataArray type=\"Int32\" Name=\"Aggregate\" format=\"ascii\">" << std::endl;
1058  fout << indent;
1059  for(int i = 0; i < int(uniqueFine.size()); i++)
1060  {
1061  fout << myAggOffset + vertex2AggId[uniqueFine[i]] << " ";
1062  if(i % 10 == 9)
1063  fout << std::endl << indent;
1064  }
1065  fout << std::endl;
1066  fout << " </DataArray>" << std::endl;
1067  fout << " <DataArray type=\"Int32\" Name=\"Processor\" format=\"ascii\">" << std::endl;
1068  fout << indent;
1069  for(int i = 0; i < int(uniqueFine.size()); i++)
1070  {
1071  fout << myRank << " ";
1072  if(i % 20 == 19)
1073  fout << std::endl << indent;
1074  }
1075  fout << std::endl;
1076  fout << " </DataArray>" << std::endl;
1077  fout << " </PointData>" << std::endl;
1078  }
1079 
1080  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1081  void VisualizationHelpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeFileVTKCoordinates(std::ofstream & fout, std::vector<int> & uniqueFine, Teuchos::ArrayRCP<const double> & fx, Teuchos::ArrayRCP<const double> & fy, Teuchos::ArrayRCP<const double> & fz, int dim) const {
1082  std::string indent = " ";
1083  fout << " <Points>" << std::endl;
1084  fout << " <DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"ascii\">" << std::endl;
1085  fout << indent;
1086  for(int i = 0; i < int(uniqueFine.size()); i++)
1087  {
1088  fout << fx[uniqueFine[i]] << " " << fy[uniqueFine[i]] << " ";
1089  if(dim == 2)
1090  fout << "0 ";
1091  else
1092  fout << fz[uniqueFine[i]] << " ";
1093  if(i % 3 == 2)
1094  fout << std::endl << indent;
1095  }
1096  fout << std::endl;
1097  fout << " </DataArray>" << std::endl;
1098  fout << " </Points>" << std::endl;
1099  }
1100 
1101  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1102  void VisualizationHelpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeFileVTKCells(std::ofstream & fout, std::vector<int> & uniqueFine, std::vector<LocalOrdinal> & vertices, std::vector<LocalOrdinal> & geomSize) const {
1103  std::string indent = " ";
1104  fout << " <Cells>" << std::endl;
1105  fout << " <DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">" << std::endl;
1106  fout << indent;
1107  for(int i = 0; i < int(vertices.size()); i++)
1108  {
1109  fout << vertices[i] << " ";
1110  if(i % 10 == 9)
1111  fout << std::endl << indent;
1112  }
1113  fout << std::endl;
1114  fout << " </DataArray>" << std::endl;
1115  fout << " <DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">" << std::endl;
1116  fout << indent;
1117  int accum = 0;
1118  for(int i = 0; i < int(geomSize.size()); i++)
1119  {
1120  accum += geomSize[i];
1121  fout << accum << " ";
1122  if(i % 10 == 9)
1123  fout << std::endl << indent;
1124  }
1125  fout << std::endl;
1126  fout << " </DataArray>" << std::endl;
1127  fout << " <DataArray type=\"Int32\" Name=\"types\" format=\"ascii\">" << std::endl;
1128  fout << indent;
1129  for(int i = 0; i < int(geomSize.size()); i++)
1130  {
1131  switch(geomSize[i])
1132  {
1133  case 1:
1134  fout << "1 "; //Point
1135  break;
1136  case 2:
1137  fout << "3 "; //Line
1138  break;
1139  case 3:
1140  fout << "5 "; //Triangle
1141  break;
1142  default:
1143  fout << "7 "; //Polygon
1144  }
1145  if(i % 30 == 29)
1146  fout << std::endl << indent;
1147  }
1148  fout << std::endl;
1149  fout << " </DataArray>" << std::endl;
1150  fout << " </Cells>" << std::endl;
1151  }
1152 
1153  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1155  fout << " </Piece>" << std::endl;
1156  fout << " </UnstructuredGrid>" << std::endl;
1157  fout << "</VTKFile>" << std::endl;
1158  }
1159 
1160 
1161  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1162  void VisualizationHelpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writePVTU(std::ofstream& pvtu, std::string baseFname, int numProcs, bool bFineEdges, bool bCoarseEdges) const {
1163  //If using vtk, filenameToWrite now contains final, correct ***.vtu filename (for the current proc)
1164  //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
1165  //pvtu file.
1166  pvtu << "<VTKFile type=\"PUnstructuredGrid\" byte_order=\"LittleEndian\">" << std::endl;
1167  pvtu << " <PUnstructuredGrid GhostLevel=\"0\">" << std::endl;
1168  pvtu << " <PPointData Scalars=\"Node Aggregate Processor\">" << std::endl;
1169  pvtu << " <PDataArray type=\"Int32\" Name=\"Node\"/>" << std::endl;
1170  pvtu << " <PDataArray type=\"Int32\" Name=\"Aggregate\"/>" << std::endl;
1171  pvtu << " <PDataArray type=\"Int32\" Name=\"Processor\"/>" << std::endl;
1172  pvtu << " </PPointData>" << std::endl;
1173  pvtu << " <PPoints>" << std::endl;
1174  pvtu << " <PDataArray type=\"Float64\" NumberOfComponents=\"3\"/>" << std::endl;
1175  pvtu << " </PPoints>" << std::endl;
1176  for(int i = 0; i < numProcs; i++) {
1177  //specify the piece for each proc (the replaceAll expression matches what the filenames will be on other procs)
1178  pvtu << " <Piece Source=\"" << replaceAll(baseFname, "%PROCID", toString(i)) << "\"/>" << std::endl;
1179  }
1180  //reference files with graph pieces, if applicable
1181  if(bFineEdges)
1182  {
1183  for(int i = 0; i < numProcs; i++)
1184  {
1185  std::string fn = replaceAll(baseFname, "%PROCID", toString(i));
1186  pvtu << " <Piece Source=\"" << fn.insert(fn.rfind(".vtu"), "-finegraph") << "\"/>" << std::endl;
1187  }
1188  }
1189  /*if(doCoarseGraphEdges_)
1190  {
1191  for(int i = 0; i < numProcs; i++)
1192  {
1193  std::string fn = replaceAll(baseFname, "%PROCID", toString(i));
1194  pvtu << " <Piece Source=\"" << fn.insert(fn.rfind(".vtu"), "-coarsegraph") << "\"/>" << std::endl;
1195  }
1196  }*/
1197  pvtu << " </PUnstructuredGrid>" << std::endl;
1198  pvtu << "</VTKFile>" << std::endl;
1199  pvtu.close();
1200  }
1201 
1202  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1204  try {
1205  std::ofstream color("random-colormap.xml");
1206  color << "<ColorMap name=\"MueLu-Random\" space=\"RGB\">" << std::endl;
1207  //Give -1, -2, -3 distinctive colors (so that the style functions can have constrasted geometry types)
1208  //Do red, orange, yellow to constrast with cool color spectrum for other types
1209  color << " <Point x=\"" << -1 << "\" o=\"1\" r=\"1\" g=\"0\" b=\"0\"/>" << std::endl;
1210  color << " <Point x=\"" << -2 << "\" o=\"1\" r=\"1\" g=\"0.6\" b=\"0\"/>" << std::endl;
1211  color << " <Point x=\"" << -3 << "\" o=\"1\" r=\"1\" g=\"1\" b=\"0\"/>" << std::endl;
1212  srand(time(NULL));
1213  for(int i = 0; i < 5000; i += 4) {
1214  color << " <Point x=\"" << i << "\" o=\"1\" r=\"" << (rand() % 50) / 256.0 << "\" g=\"" << (rand() % 256) / 256.0 << "\" b=\"" << (rand() % 256) / 256.0 << "\"/>" << std::endl;
1215  }
1216  color << "</ColorMap>" << std::endl;
1217  color.close();
1218  }
1219  catch(std::exception& e) {
1220  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError,
1221  "VisualizationHelpers::buildColormap: Error while building colormap file: " << e.what());
1222  }
1223  }
1224 
1225 } // namespace MueLu
1226 
1227 #endif /* MUELU_VISUALIZATIONHELPERS_DEF_HPP_ */
static void doPointCloud(std::vector< int > &vertices, std::vector< int > &geomSizes, LO numLocalAggs, LO numFineNodes)
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
std::string getPVTUFileName(int numProcs, int myRank, int level, const Teuchos::ParameterList &pL) const
void writeFileVTKOpening(std::ofstream &fout, std::vector< int > &uniqueFine, std::vector< int > &geomSizesFine) const
std::string replaceAll(std::string result, const std::string &replaceWhat, const std::string &replaceWithWhat) const
void writeFileVTKData(std::ofstream &fout, std::vector< int > &uniqueFine, LocalOrdinal myAggOffset, ArrayRCP< LocalOrdinal > &vertex2AggId, int myRank) const
std::vector< int > makeUnique(std::vector< int > &vertices) const
replaces node indices in vertices with compressed unique indices, and returns list of unique points ...
std::string getFileName(int numProcs, int myRank, int level, const Teuchos::ParameterList &pL) const
Namespace for MueLu classes and methods.
void writeFileVTKNodes(std::ofstream &fout, std::vector< int > &uniqueFine, Teuchos::RCP< const Map > &nodeMap) const
RCP< ParameterList > GetValidParameterList() const
static double pointDistFromTri(myVec3 point, myVec3 v1, myVec3 v2, myVec3 v3)
static double dotProduct(myVec2 v1, myVec2 v2)
void writePVTU(std::ofstream &pvtu, std::string baseFname, int numProcs, bool bFineEdges=false, bool bCoarseEdges=false) const
std::string getBaseFileName(int numProcs, int level, const Teuchos::ParameterList &pL) const
static void doGraphEdges(std::vector< int > &vertices, std::vector< int > &geomSizes, Teuchos::RCP< GraphBase > &G, Teuchos::ArrayRCP< const double > &fx, Teuchos::ArrayRCP< const double > &fy, Teuchos::ArrayRCP< const double > &fz)
static std::vector< int > giftWrap(std::vector< myVec2 > &points, std::vector< int > &nodes, const Teuchos::ArrayRCP< const double > &xCoords, const Teuchos::ArrayRCP< const double > &yCoords)
static void doJacks(std::vector< int > &vertices, std::vector< int > &geomSizes, LO numLocalAggs, LO numFineNodes, const std::vector< bool > &isRoot, const ArrayRCP< LO > &vertex2AggId)
static void doConvexHulls2D(std::vector< int > &vertices, std::vector< int > &geomSizes, LO numLocalAggs, LO numFineNodes, const std::vector< bool > &isRoot, const ArrayRCP< LO > &vertex2AggId, const Teuchos::ArrayRCP< const double > &xCoords, const Teuchos::ArrayRCP< const double > &yCoords, const Teuchos::ArrayRCP< const double > &zCoords)
static myVec3 crossProduct(myVec3 v1, myVec3 v2)
void writeFileVTKCoordinates(std::ofstream &fout, std::vector< int > &uniqueFine, Teuchos::ArrayRCP< const double > &fx, Teuchos::ArrayRCP< const double > &fy, Teuchos::ArrayRCP< const double > &fz, int dim) const
static myVec2 vecSubtract(myVec2 v1, myVec2 v2)
static std::vector< myTriangle > processTriangle(std::list< myTriangle > &tris, myTriangle tri, std::list< int > &pointsInFront, myVec3 &barycenter, const Teuchos::ArrayRCP< const double > &xCoords, const Teuchos::ArrayRCP< const double > &yCoords, const Teuchos::ArrayRCP< const double > &zCoords)
static void doConvexHulls3D(std::vector< int > &vertices, std::vector< int > &geomSizes, LO numLocalAggs, LO numFineNodes, const std::vector< bool > &isRoot, const ArrayRCP< LO > &vertex2AggId, const Teuchos::ArrayRCP< const double > &xCoords, const Teuchos::ArrayRCP< const double > &yCoords, const Teuchos::ArrayRCP< const double > &zCoords)
static double distance(myVec2 p1, myVec2 p2)
void writeFileVTKCells(std::ofstream &fout, std::vector< int > &uniqueFine, std::vector< LocalOrdinal > &vertices, std::vector< LocalOrdinal > &geomSize) const
static bool isInFront(myVec3 point, myVec3 inPlane, myVec3 n)
Exception throws to report errors in the internal logical of the program.
void replaceAll(std::string &str, const std::string &from, const std::string &to)
void writeFileVTKClosing(std::ofstream &fout) const