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 #ifdef HAVE_MUELU_CGAL
64 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
65 #include <CGAL/convex_hull_2.h>
66 #include <CGAL/Polyhedron_3.h>
67 #include <CGAL/convex_hull_3.h>
68 #endif
69 
70 
71 namespace MueLu {
72 
73  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
75  RCP<ParameterList> validParamList = rcp(new ParameterList());
76 
77  validParamList->set< std::string > ("visualization: output filename", "viz%LEVELID", "filename for VTK-style visualization output");
78  validParamList->set< int > ("visualization: output file: time step", 0, "time step variable for output file name");// Remove me?
79  validParamList->set< int > ("visualization: output file: iter", 0, "nonlinear iteration variable for output file name");//Remove me?
80  validParamList->set<std::string> ("visualization: style", "Point Cloud", "style of aggregate visualization for VTK output. Can be 'Point Cloud', 'Jacks', 'Convex Hulls'");
81  validParamList->set<bool> ("visualization: build colormap", false, "Whether to build a random color map in a separate xml file.");
82  validParamList->set<bool> ("visualization: fine graph edges", false, "Whether to draw all fine node connections along with the aggregates.");
83 
84  return validParamList;
85  }
86 
87  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
88  void VisualizationHelpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::doPointCloud(std::vector<int>& vertices, std::vector<int>& geomSizes, LO /* numLocalAggs */, LO numFineNodes) {
89  vertices.reserve(numFineNodes);
90  geomSizes.reserve(numFineNodes);
91  for(LocalOrdinal i = 0; i < numFineNodes; i++)
92  {
93  vertices.push_back(i);
94  geomSizes.push_back(1);
95  }
96  }
97 
98  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
99  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) {
100  //For each aggregate, find the root node then connect it with the other nodes in the aggregate
101  //Doesn't matter the order, as long as all edges are found.
102  vertices.reserve(vertices.size() + 3 * (numFineNodes - numLocalAggs));
103  geomSizes.reserve(vertices.size() + 2 * (numFineNodes - numLocalAggs));
104  int root = 0;
105  for(int i = 0; i < numLocalAggs; i++) //TODO: Replace this O(n^2) with a better way
106  {
107  while(!isRoot[root])
108  root++;
109  int numInAggFound = 0;
110  for(int j = 0; j < numFineNodes; j++)
111  {
112  if(j == root) //don't make a connection from the root to itself
113  {
114  numInAggFound++;
115  continue;
116  }
117  if(vertex2AggId[root] == vertex2AggId[j])
118  {
119  vertices.push_back(root);
120  vertices.push_back(j);
121  geomSizes.push_back(2);
122  //Also draw the free endpoint explicitly for the current line
123  vertices.push_back(j);
124  geomSizes.push_back(1);
125  numInAggFound++;
126  //if(numInAggFound == aggSizes_[vertex2AggId_[root]]) //don't spend more time looking if done with that root
127  // break;
128  }
129  }
130  root++; //get set up to look for the next root
131  }
132  }
133 
134  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
135  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 typename Teuchos::ScalarTraits<Scalar>::coordinateType>& xCoords, const Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType>& yCoords) {
136 
137  // This algorithm is based on Andrew's Monotone Chain variant of the Graham Scan for Convex Hulls. It adds
138  // a colinearity check which we'll need for our viz.
139  for(int agg = 0; agg < numLocalAggs; agg++) {
140  std::vector<int> aggNodes;
141  for(int i = 0; i < numFineNodes; i++) {
142  if(vertex2AggId[i] == agg)
143  aggNodes.push_back(i);
144  }
145 
146  //have a list of nodes in the aggregate
147  TEUCHOS_TEST_FOR_EXCEPTION(aggNodes.size() == 0, Exceptions::RuntimeError,
148  "CoarseningVisualization::doConvexHulls2D: aggregate contains zero nodes!");
149  if(aggNodes.size() == 1) {
150  vertices.push_back(aggNodes.front());
151  geomSizes.push_back(1);
152  continue;
153  }
154  if(aggNodes.size() == 2) {
155  vertices.push_back(aggNodes.front());
156  vertices.push_back(aggNodes.back());
157  geomSizes.push_back(2);
158  continue;
159  }
160  else {
161  int N = (int) aggNodes.size();
162  using MyPair = std::pair<myVec2,int>;
163  std::vector<MyPair> pointsAndIndex(N);
164  for(int i=0; i<N; i++) {
165  pointsAndIndex[i] = std::make_pair(myVec2(xCoords[aggNodes[i]], yCoords[aggNodes[i]]),i);
166  }
167 
168  // Sort by x coordinate
169  std::sort(pointsAndIndex.begin(),pointsAndIndex.end(),[](const MyPair &a, const MyPair &b) {
170  return a.first.x < b.first.x || (a.first.x == b.first.x && a.first.y < b.first.y);
171  });
172 
173 
174  // Colinearity check
175  bool colinear=true;
176  for(int i=0; i<N; i++) {
177  if(ccw(pointsAndIndex[i].first,pointsAndIndex[(i+1)%N].first,pointsAndIndex[(i+2)%N].first)!=0) {
178  colinear=false;
179  break;
180  }
181  }
182 
183  if(colinear) {
184  vertices.push_back(aggNodes[pointsAndIndex.front().second]);
185  vertices.push_back(aggNodes[pointsAndIndex.back().second]);
186  geomSizes.push_back(2);
187  }
188  else {
189  std::vector<int> hull(2*N);
190  int count=0;
191 
192  // Build lower hull
193  for(int i=0; i<N; i++) {
194  while (count >=2 && ccw(pointsAndIndex[hull[count-2]].first,pointsAndIndex[hull[count-1]].first,pointsAndIndex[i].first)<=0) {
195  count--;
196  }
197  hull[count++] = i;
198  }
199 
200  // Build the upper hull
201  int t=count+1;
202  for(int i=N-1; i>0; i--) {
203  while(count >= t && ccw(pointsAndIndex[hull[count-2]].first,pointsAndIndex[hull[count-1]].first,pointsAndIndex[i-1].first)<=0) {
204  count--;
205  }
206  hull[count++]=i-1;
207  }
208 
209  // Remove the duplicated point
210  hull.resize(count-1);
211 
212  // Verify: Check that hull retains CCW order
213  for(int i=0; i<(int)hull.size(); i++) {
214  TEUCHOS_TEST_FOR_EXCEPTION(ccw(pointsAndIndex[hull[i]].first,pointsAndIndex[hull[(i+1)%hull.size()]].first,pointsAndIndex[hull[(i+1)%hull.size()]].first) == 1, Exceptions::RuntimeError,"CoarseningVisualization::doConvexHulls2D: counter-clockwise verification fails");
215  }
216 
217  // We now need to undo the indices from the sorting
218  for(int i=0; i<(int)hull.size(); i++) {
219  hull[i] = aggNodes[pointsAndIndex[hull[i]].second];
220  }
221 
222  vertices.reserve(vertices.size() + hull.size());
223  for(size_t i = 0; i < hull.size(); i++) {
224  vertices.push_back(hull[i]);
225  }
226  geomSizes.push_back(hull.size());
227  }//else colinear
228  }//else 3 + nodes
229  }
230  }
231 
232 #ifdef HAVE_MUELU_CGAL
233  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
234  void VisualizationHelpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::doCGALConvexHulls2D(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 typename Teuchos::ScalarTraits<Scalar>::coordinateType>& xCoords, const Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType>& yCoords) {
235 
236  typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
237  typedef K::Point_2 Point_2;
238 
239  for(int agg = 0; agg < numLocalAggs; agg++) {
240  std::vector<int> aggNodes;
241  std::vector<Point_2> aggPoints;
242  for(int i = 0; i < numFineNodes; i++) {
243  if(vertex2AggId[i] == agg) {
244  Point_2 p(xCoords[i], yCoords[i]);
245  aggPoints.push_back(p);
246  aggNodes.push_back(i);
247  }
248  }
249  //have a list of nodes in the aggregate
250  TEUCHOS_TEST_FOR_EXCEPTION(aggNodes.size() == 0, Exceptions::RuntimeError,
251  "CoarseningVisualization::doCGALConvexHulls2D: aggregate contains zero nodes!");
252  if(aggNodes.size() == 1) {
253  vertices.push_back(aggNodes.front());
254  geomSizes.push_back(1);
255  continue;
256  }
257  if(aggNodes.size() == 2) {
258  vertices.push_back(aggNodes.front());
259  vertices.push_back(aggNodes.back());
260  geomSizes.push_back(2);
261  continue;
262  }
263  //check if all points are collinear, need to explicitly draw a line in that case.
264  bool collinear = true; //assume true at first, if a segment not parallel to others then clear
265  {
266  std::vector<int>::iterator it = aggNodes.begin();
267  myVec3 firstPoint(xCoords[*it], yCoords[*it], 0);
268  it++;
269  myVec3 secondPoint(xCoords[*it], yCoords[*it], 0);
270  it++; //it now points to third node in the aggregate
271  myVec3 norm1(-(secondPoint.y - firstPoint.y), secondPoint.x - firstPoint.x, 0);
272  do {
273  myVec3 thisNorm(yCoords[*it] - firstPoint.y, firstPoint.x - xCoords[*it], 0);
274  //rotate one of the vectors by 90 degrees so that dot product is 0 if the two are parallel
275  double temp = thisNorm.x;
276  thisNorm.x = thisNorm.y;
277  thisNorm.y = temp;
278  double comp = dotProduct(norm1, thisNorm);
279  if(-1e-8 > comp || comp > 1e-8) {
280  collinear = false;
281  break;
282  }
283  it++;
284  }
285  while(it != aggNodes.end());
286  }
287  if(collinear)
288  {
289  //find the most distant two points in the plane and use as endpoints of line representing agg
290  std::vector<int>::iterator min = aggNodes.begin(); //min X then min Y where x is a tie
291  std::vector<int>::iterator max = aggNodes.begin(); //max X then max Y where x is a tie
292  for(std::vector<int>::iterator it = ++aggNodes.begin(); it != aggNodes.end(); it++) {
293  if(xCoords[*it] < xCoords[*min])
294  min = it;
295  else if(xCoords[*it] == xCoords[*min]) {
296  if(yCoords[*it] < yCoords[*min])
297  min = it;
298  }
299  if(xCoords[*it] > xCoords[*max])
300  max = it;
301  else if(xCoords[*it] == xCoords[*max]) {
302  if(yCoords[*it] > yCoords[*max])
303  max = it;
304  }
305  }
306  //Just set up a line between nodes *min and *max
307  vertices.push_back(*min);
308  vertices.push_back(*max);
309  geomSizes.push_back(2);
310  continue; //jump to next aggregate in loop
311  }
312  // aggregate has more than 2 points and is not collinear
313  {
314  std::vector<Point_2> result;
315  CGAL::convex_hull_2( aggPoints.begin(), aggPoints.end(), std::back_inserter(result) );
316  // loop over all result points and find the corresponding node id
317  for (size_t r = 0; r < result.size(); r++) {
318  // loop over all aggregate nodes and find corresponding node id
319  for(size_t l = 0; l < aggNodes.size(); l++)
320  {
321  if(fabs(result[r].x() - xCoords[aggNodes[l]]) < 1e-12 &&
322  fabs(result[r].y() - yCoords[aggNodes[l]]) < 1e-12)
323  {
324  vertices.push_back(aggNodes[l]);
325  break;
326  }
327  }
328 
329  }
330  geomSizes.push_back(result.size());
331  }
332  }
333  }
334 
335 #endif
336 
337  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
338  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 typename Teuchos::ScalarTraits<Scalar>::coordinateType>& xCoords, const Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType>& yCoords, const Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType>& zCoords) {
339  //Use 3D quickhull algo.
340  //Vector of node indices representing triangle vertices
341  //Note: Calculate the hulls first since will only include point data for points in the hulls
342  //Effectively the size() of vertIndices after each hull is added to it
343  typedef std::list<int>::iterator Iter;
344  for(int agg = 0; agg < numLocalAggs; agg++) {
345  std::list<int> aggNodes; //At first, list of all nodes in the aggregate. As nodes are enclosed or included by/in hull, remove them
346  for(int i = 0; i < numFineNodes; i++) {
347  if(vertex2AggId[i] == agg)
348  aggNodes.push_back(i);
349  }
350  //First, check anomalous cases
351  TEUCHOS_TEST_FOR_EXCEPTION(aggNodes.size() == 0, Exceptions::RuntimeError,
352  "CoarseningVisualization::doConvexHulls3D: aggregate contains zero nodes!");
353  if(aggNodes.size() == 1) {
354  vertices.push_back(aggNodes.front());
355  geomSizes.push_back(1);
356  continue;
357  } else if(aggNodes.size() == 2) {
358  vertices.push_back(aggNodes.front());
359  vertices.push_back(aggNodes.back());
360  geomSizes.push_back(2);
361  continue;
362  }
363  //check for collinearity
364  bool areCollinear = true;
365  {
366  Iter it = aggNodes.begin();
367  myVec3 firstVec(xCoords[*it], yCoords[*it], zCoords[*it]);
368  myVec3 comp;
369  {
370  it++;
371  myVec3 secondVec(xCoords[*it], yCoords[*it], zCoords[*it]); //cross this with other vectors to compare
372  comp = vecSubtract(secondVec, firstVec);
373  it++;
374  }
375  for(; it != aggNodes.end(); it++) {
376  myVec3 thisVec(xCoords[*it], yCoords[*it], zCoords[*it]);
377  myVec3 cross = crossProduct(vecSubtract(thisVec, firstVec), comp);
378  if(mymagnitude(cross) > 1e-10) {
379  areCollinear = false;
380  break;
381  }
382  }
383  }
384  if(areCollinear)
385  {
386  //find the endpoints of segment describing all the points
387  //compare x, if tie compare y, if tie compare z
388  Iter min = aggNodes.begin();
389  Iter max = aggNodes.begin();
390  Iter it = ++aggNodes.begin();
391  for(; it != aggNodes.end(); it++) {
392  if(xCoords[*it] < xCoords[*min]) min = it;
393  else if(xCoords[*it] == xCoords[*min]) {
394  if(yCoords[*it] < yCoords[*min]) min = it;
395  else if(yCoords[*it] == yCoords[*min]) {
396  if(zCoords[*it] < zCoords[*min]) min = it;
397  }
398  }
399  if(xCoords[*it] > xCoords[*max]) max = it;
400  else if(xCoords[*it] == xCoords[*max]) {
401  if(yCoords[*it] > yCoords[*max]) max = it;
402  else if(yCoords[*it] == yCoords[*max]) {
403  if(zCoords[*it] > zCoords[*max])
404  max = it;
405  }
406  }
407  }
408  vertices.push_back(*min);
409  vertices.push_back(*max);
410  geomSizes.push_back(2);
411  continue;
412  }
413  bool areCoplanar = true;
414  {
415  //number of points is known to be >= 3
416  Iter vert = aggNodes.begin();
417  myVec3 v1(xCoords[*vert], yCoords[*vert], zCoords[*vert]);
418  vert++;
419  myVec3 v2(xCoords[*vert], yCoords[*vert], zCoords[*vert]);
420  vert++;
421  myVec3 v3(xCoords[*vert], yCoords[*vert], zCoords[*vert]);
422  vert++;
423  //Make sure the first three points aren't also collinear (need a non-degenerate triangle to get a normal)
424  while(mymagnitude(crossProduct(vecSubtract(v1, v2), vecSubtract(v1, v3))) < 1e-10) {
425  //Replace the third point with the next point
426  v3 = myVec3(xCoords[*vert], yCoords[*vert], zCoords[*vert]);
427  vert++;
428  }
429  for(; vert != aggNodes.end(); vert++) {
430  myVec3 pt(xCoords[*vert], yCoords[*vert], zCoords[*vert]);
431  if(fabs(pointDistFromTri(pt, v1, v2, v3)) > 1e-12) {
432  areCoplanar = false;
433  break;
434  }
435  }
436  if(areCoplanar) {
437  //do 2D convex hull
438  myVec3 planeNorm = getNorm(v1, v2, v3);
439  planeNorm.x = fabs(planeNorm.x);
440  planeNorm.y = fabs(planeNorm.y);
441  planeNorm.z = fabs(planeNorm.z);
442  std::vector<myVec2> points;
443  std::vector<int> nodes;
444  if(planeNorm.x >= planeNorm.y && planeNorm.x >= planeNorm.z) {
445  //project points to yz plane to make hull
446  for(Iter it = aggNodes.begin(); it != aggNodes.end(); it++) {
447  nodes.push_back(*it);
448  points.push_back(myVec2(yCoords[*it], zCoords[*it]));
449  }
450  }
451  if(planeNorm.y >= planeNorm.x && planeNorm.y >= planeNorm.z) {
452  //use xz
453  for(Iter it = aggNodes.begin(); it != aggNodes.end(); it++) {
454  nodes.push_back(*it);
455  points.push_back(myVec2(xCoords[*it], zCoords[*it]));
456  }
457  }
458  if(planeNorm.z >= planeNorm.x && planeNorm.z >= planeNorm.y) {
459  for(Iter it = aggNodes.begin(); it != aggNodes.end(); it++) {
460  nodes.push_back(*it);
461  points.push_back(myVec2(xCoords[*it], yCoords[*it]));
462  }
463  }
464  std::vector<int> convhull2d = giftWrap(points, nodes, xCoords, yCoords);
465  geomSizes.push_back(convhull2d.size());
466  vertices.reserve(vertices.size() + convhull2d.size());
467  for(size_t i = 0; i < convhull2d.size(); i++)
468  vertices.push_back(convhull2d[i]);
469  continue;
470  }
471  }
472  Iter exIt = aggNodes.begin(); //iterator to be used for searching for min/max x/y/z
473  int extremeSix[] = {*exIt, *exIt, *exIt, *exIt, *exIt, *exIt}; //nodes with minimumX, maxX, minY ...
474  exIt++;
475  for(; exIt != aggNodes.end(); exIt++) {
476  Iter it = exIt;
477  if(xCoords[*it] < xCoords[extremeSix[0]] ||
478  (xCoords[*it] == xCoords[extremeSix[0]] && yCoords[*it] < yCoords[extremeSix[0]]) ||
479  (xCoords[*it] == xCoords[extremeSix[0]] && yCoords[*it] == yCoords[extremeSix[0]] && zCoords[*it] < zCoords[extremeSix[0]]))
480  extremeSix[0] = *it;
481  if(xCoords[*it] > xCoords[extremeSix[1]] ||
482  (xCoords[*it] == xCoords[extremeSix[1]] && yCoords[*it] > yCoords[extremeSix[1]]) ||
483  (xCoords[*it] == xCoords[extremeSix[1]] && yCoords[*it] == yCoords[extremeSix[1]] && zCoords[*it] > zCoords[extremeSix[1]]))
484  extremeSix[1] = *it;
485  if(yCoords[*it] < yCoords[extremeSix[2]] ||
486  (yCoords[*it] == yCoords[extremeSix[2]] && zCoords[*it] < zCoords[extremeSix[2]]) ||
487  (yCoords[*it] == yCoords[extremeSix[2]] && zCoords[*it] == zCoords[extremeSix[2]] && xCoords[*it] < xCoords[extremeSix[2]]))
488  extremeSix[2] = *it;
489  if(yCoords[*it] > yCoords[extremeSix[3]] ||
490  (yCoords[*it] == yCoords[extremeSix[3]] && zCoords[*it] > zCoords[extremeSix[3]]) ||
491  (yCoords[*it] == yCoords[extremeSix[3]] && zCoords[*it] == zCoords[extremeSix[3]] && xCoords[*it] > xCoords[extremeSix[3]]))
492  extremeSix[3] = *it;
493  if(zCoords[*it] < zCoords[extremeSix[4]] ||
494  (zCoords[*it] == zCoords[extremeSix[4]] && xCoords[*it] < xCoords[extremeSix[4]]) ||
495  (zCoords[*it] == zCoords[extremeSix[4]] && xCoords[*it] == xCoords[extremeSix[4]] && yCoords[*it] < yCoords[extremeSix[4]]))
496  extremeSix[4] = *it;
497  if(zCoords[*it] > zCoords[extremeSix[5]] ||
498  (zCoords[*it] == zCoords[extremeSix[5]] && xCoords[*it] > xCoords[extremeSix[5]]) ||
499  (zCoords[*it] == zCoords[extremeSix[5]] && xCoords[*it] == xCoords[extremeSix[5]] && yCoords[*it] > zCoords[extremeSix[5]]))
500  extremeSix[5] = *it;
501  }
502  myVec3 extremeVectors[6];
503  for(int i = 0; i < 6; i++) {
504  myVec3 thisExtremeVec(xCoords[extremeSix[i]], yCoords[extremeSix[i]], zCoords[extremeSix[i]]);
505  extremeVectors[i] = thisExtremeVec;
506  }
507  double maxDist = 0;
508  int base1 = 0; //ints from 0-5: which pair out of the 6 extreme points are the most distant? (indices in extremeSix and extremeVectors)
509  int base2 = 0;
510  for(int i = 0; i < 5; i++) {
511  for(int j = i + 1; j < 6; j++) {
512  double thisDist = distance(extremeVectors[i], extremeVectors[j]);
513  // Want to make sure thisDist is significantly larger than maxDist to prevent roundoff errors from impacting node choice.
514  if(thisDist > maxDist && thisDist - maxDist > 1e-12) {
515  maxDist = thisDist;
516  base1 = i;
517  base2 = j;
518  }
519  }
520  }
521  std::list<myTriangle> hullBuilding; //each Triangle is a triplet of nodes (int IDs) that form a triangle
522  //remove base1 and base2 iters from aggNodes, they are known to be in the hull
523  aggNodes.remove(extremeSix[base1]);
524  aggNodes.remove(extremeSix[base2]);
525  //extremeSix[base1] and [base2] still have the myVec3 representation
526  myTriangle tri;
527  tri.v1 = extremeSix[base1];
528  tri.v2 = extremeSix[base2];
529  //Now find the point that is furthest away from the line between base1 and base2
530  maxDist = 0;
531  //need the vectors to do "quadruple product" formula
532  myVec3 b1 = extremeVectors[base1];
533  myVec3 b2 = extremeVectors[base2];
534  Iter thirdNode;
535  for(Iter node = aggNodes.begin(); node != aggNodes.end(); node++) {
536  myVec3 nodePos(xCoords[*node], yCoords[*node], zCoords[*node]);
537  double dist = mymagnitude(crossProduct(vecSubtract(nodePos, b1), vecSubtract(nodePos, b2))) / mymagnitude(vecSubtract(b2, b1));
538  // Want to make sure dist is significantly larger than maxDist to prevent roundoff errors from impacting node choice.
539  if(dist > maxDist && dist - maxDist > 1e-12) {
540  maxDist = dist;
541  thirdNode = node;
542  }
543  }
544  //Now know the last node in the first triangle
545  tri.v3 = *thirdNode;
546  hullBuilding.push_back(tri);
547  myVec3 b3(xCoords[*thirdNode], yCoords[*thirdNode], zCoords[*thirdNode]);
548  aggNodes.erase(thirdNode);
549  //Find the fourth node (most distant from triangle) to form tetrahedron
550  maxDist = 0;
551  int fourthVertex = -1;
552  for(Iter node = aggNodes.begin(); node != aggNodes.end(); node++) {
553  myVec3 thisNode(xCoords[*node], yCoords[*node], zCoords[*node]);
554  double nodeDist = pointDistFromTri(thisNode, b1, b2, b3);
555  // Want to make sure nodeDist is significantly larger than maxDist to prevent roundoff errors from impacting node choice.
556  if(nodeDist > maxDist && nodeDist - maxDist > 1e-12) {
557  maxDist = nodeDist;
558  fourthVertex = *node;
559  }
560  }
561  aggNodes.remove(fourthVertex);
562  myVec3 b4(xCoords[fourthVertex], yCoords[fourthVertex], zCoords[fourthVertex]);
563  //Add three new triangles to hullBuilding to form the first tetrahedron
564  //use tri to hold the triangle info temporarily before being added to list
565  tri = hullBuilding.front();
566  tri.v1 = fourthVertex;
567  hullBuilding.push_back(tri);
568  tri = hullBuilding.front();
569  tri.v2 = fourthVertex;
570  hullBuilding.push_back(tri);
571  tri = hullBuilding.front();
572  tri.v3 = fourthVertex;
573  hullBuilding.push_back(tri);
574  //now orient all four triangles so that the vertices are oriented clockwise (so getNorm_ points outward for each)
575  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);
576  for(std::list<myTriangle>::iterator tetTri = hullBuilding.begin(); tetTri != hullBuilding.end(); tetTri++) {
577  myVec3 triVert1(xCoords[tetTri->v1], yCoords[tetTri->v1], zCoords[tetTri->v1]);
578  myVec3 triVert2(xCoords[tetTri->v2], yCoords[tetTri->v2], zCoords[tetTri->v2]);
579  myVec3 triVert3(xCoords[tetTri->v3], yCoords[tetTri->v3], zCoords[tetTri->v3]);
580  myVec3 trinorm = getNorm(triVert1, triVert2, triVert3);
581  myVec3 ptInPlane = tetTri == hullBuilding.begin() ? b1 : b4; //first triangle definitely has b1 in it, other three definitely have b4
582  if(isInFront(barycenter, ptInPlane, trinorm)) {
583  //don't want the faces of the tetrahedron to face inwards (towards barycenter) so reverse orientation
584  //by swapping two vertices
585  int temp = tetTri->v1;
586  tetTri->v1 = tetTri->v2;
587  tetTri->v2 = temp;
588  }
589  }
590  //now, have starting polyhedron in hullBuilding (all faces are facing outwards according to getNorm_) and remaining nodes to process are in aggNodes
591  //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.
592  //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
593  //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
594  //new triangles, since they are now enclosed in the hull.
595  //Construct point lists for each face of the tetrahedron individually.
596  myVec3 trinorms[4]; //normals to the four tetrahedron faces, now oriented outwards
597  int index = 0;
598  for(std::list<myTriangle>::iterator tetTri = hullBuilding.begin(); tetTri != hullBuilding.end(); tetTri++) {
599  myVec3 triVert1(xCoords[tetTri->v1], yCoords[tetTri->v1], zCoords[tetTri->v1]);
600  myVec3 triVert2(xCoords[tetTri->v2], yCoords[tetTri->v2], zCoords[tetTri->v2]);
601  myVec3 triVert3(xCoords[tetTri->v3], yCoords[tetTri->v3], zCoords[tetTri->v3]);
602  trinorms[index] = getNorm(triVert1, triVert2, triVert3);
603  index++;
604  }
605  std::list<int> startPoints1;
606  std::list<int> startPoints2;
607  std::list<int> startPoints3;
608  std::list<int> startPoints4;
609  //scope this so that 'point' is not in function scope
610  {
611  Iter point = aggNodes.begin();
612  while(!aggNodes.empty()) //this removes points one at a time as they are put in startPointsN or are already done
613  {
614  myVec3 pointVec(xCoords[*point], yCoords[*point], zCoords[*point]);
615  //Note: Because of the way the tetrahedron faces are constructed above,
616  //face 1 definitely contains b1 and faces 2-4 definitely contain b4.
617  if(isInFront(pointVec, b1, trinorms[0])) {
618  startPoints1.push_back(*point);
619  point = aggNodes.erase(point);
620  } else if(isInFront(pointVec, b4, trinorms[1])) {
621  startPoints2.push_back(*point);
622  point = aggNodes.erase(point);
623  } else if(isInFront(pointVec, b4, trinorms[2])) {
624  startPoints3.push_back(*point);
625  point = aggNodes.erase(point);
626  } else if(isInFront(pointVec, b4, trinorms[3])) {
627  startPoints4.push_back(*point);
628  point = aggNodes.erase(point);
629  } else {
630  point = aggNodes.erase(point); //points here are already inside tetrahedron.
631  }
632  }
633  //Call processTriangle for each triangle in the initial tetrahedron, one at a time.
634  }
635  typedef std::list<myTriangle>::iterator TriIter;
636  TriIter firstTri = hullBuilding.begin();
637  myTriangle start1 = *firstTri;
638  firstTri++;
639  myTriangle start2 = *firstTri;
640  firstTri++;
641  myTriangle start3 = *firstTri;
642  firstTri++;
643  myTriangle start4 = *firstTri;
644  //kick off depth-first recursive filling of hullBuilding list with all triangles in the convex hull
645  if(!startPoints1.empty())
646  processTriangle(hullBuilding, start1, startPoints1, barycenter, xCoords, yCoords, zCoords);
647  if(!startPoints2.empty())
648  processTriangle(hullBuilding, start2, startPoints2, barycenter, xCoords, yCoords, zCoords);
649  if(!startPoints3.empty())
650  processTriangle(hullBuilding, start3, startPoints3, barycenter, xCoords, yCoords, zCoords);
651  if(!startPoints4.empty())
652  processTriangle(hullBuilding, start4, startPoints4, barycenter, xCoords, yCoords, zCoords);
653  //hullBuilding now has all triangles that make up this hull.
654  //Dump hullBuilding info into the list of all triangles for the scene.
655  vertices.reserve(vertices.size() + 3 * hullBuilding.size());
656  for(TriIter hullTri = hullBuilding.begin(); hullTri != hullBuilding.end(); hullTri++) {
657  vertices.push_back(hullTri->v1);
658  vertices.push_back(hullTri->v2);
659  vertices.push_back(hullTri->v3);
660  geomSizes.push_back(3);
661  }
662  }
663  }
664 
665 #ifdef HAVE_MUELU_CGAL
666  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
667  void VisualizationHelpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::doCGALConvexHulls3D(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 typename Teuchos::ScalarTraits<Scalar>::coordinateType>& xCoords, const Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType>& yCoords, const Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType>& zCoords) {
668 
669  typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
670  typedef K::Point_3 Point_3;
671  typedef CGAL::Polyhedron_3<K> Polyhedron_3;
672  typedef std::vector<int>::iterator Iter;
673  for(int agg = 0; agg < numLocalAggs; agg++) {
674  std::vector<int> aggNodes;
675  std::vector<Point_3> aggPoints;
676  for(int i = 0; i < numFineNodes; i++) {
677  if(vertex2AggId[i] == agg) {
678  Point_3 p(xCoords[i], yCoords[i], zCoords[i]);
679  aggPoints.push_back(p);
680  aggNodes.push_back(i);
681  }
682  }
683  //First, check anomalous cases
684  TEUCHOS_TEST_FOR_EXCEPTION(aggNodes.size() == 0, Exceptions::RuntimeError,
685  "CoarseningVisualization::doCGALConvexHulls3D: aggregate contains zero nodes!");
686  if(aggNodes.size() == 1) {
687  vertices.push_back(aggNodes.front());
688  geomSizes.push_back(1);
689  continue;
690  } else if(aggNodes.size() == 2) {
691  vertices.push_back(aggNodes.front());
692  vertices.push_back(aggNodes.back());
693  geomSizes.push_back(2);
694  continue;
695  }
696  //check for collinearity
697  bool areCollinear = true;
698  {
699  Iter it = aggNodes.begin();
700  myVec3 firstVec(xCoords[*it], yCoords[*it], zCoords[*it]);
701  myVec3 comp;
702  {
703  it++;
704  myVec3 secondVec(xCoords[*it], yCoords[*it], zCoords[*it]); //cross this with other vectors to compare
705  comp = vecSubtract(secondVec, firstVec);
706  it++;
707  }
708  for(; it != aggNodes.end(); it++) {
709  myVec3 thisVec(xCoords[*it], yCoords[*it], zCoords[*it]);
710  myVec3 cross = crossProduct(vecSubtract(thisVec, firstVec), comp);
711  if(mymagnitude(cross) > 1e-8) {
712  areCollinear = false;
713  break;
714  }
715  }
716  }
717  if(areCollinear)
718  {
719  //find the endpoints of segment describing all the points
720  //compare x, if tie compare y, if tie compare z
721  Iter min = aggNodes.begin();
722  Iter max = aggNodes.begin();
723  Iter it = ++aggNodes.begin();
724  for(; it != aggNodes.end(); it++) {
725  if(xCoords[*it] < xCoords[*min]) min = it;
726  else if(xCoords[*it] == xCoords[*min]) {
727  if(yCoords[*it] < yCoords[*min]) min = it;
728  else if(yCoords[*it] == yCoords[*min]) {
729  if(zCoords[*it] < zCoords[*min]) min = it;
730  }
731  }
732  if(xCoords[*it] > xCoords[*max]) max = it;
733  else if(xCoords[*it] == xCoords[*max]) {
734  if(yCoords[*it] > yCoords[*max]) max = it;
735  else if(yCoords[*it] == yCoords[*max]) {
736  if(zCoords[*it] > zCoords[*max])
737  max = it;
738  }
739  }
740  }
741  vertices.push_back(*min);
742  vertices.push_back(*max);
743  geomSizes.push_back(2);
744  continue;
745  }
746  // do not handle coplanar or general case here. Just let's use CGAL
747  {
748  Polyhedron_3 result;
749  CGAL::convex_hull_3( aggPoints.begin(), aggPoints.end(), result);
750 
751  // loop over all facets
752  Polyhedron_3::Facet_iterator fi;
753  for (fi = result.facets_begin(); fi != result.facets_end(); fi++) {
754  int cntVertInAgg = 0;
755  Polyhedron_3::Halfedge_around_facet_const_circulator hit = fi->facet_begin();
756  do {
757  const Point_3 & pp = hit->vertex()->point();
758  // loop over all aggregate nodes and find corresponding node id
759  for(size_t l = 0; l < aggNodes.size(); l++)
760  {
761  if(fabs(pp.x() - xCoords[aggNodes[l]]) < 1e-12 &&
762  fabs(pp.y() - yCoords[aggNodes[l]]) < 1e-12 &&
763  fabs(pp.z() - zCoords[aggNodes[l]]) < 1e-12)
764  {
765  vertices.push_back(aggNodes[l]);
766  cntVertInAgg++;
767  break;
768  }
769  }
770  } while (++hit != fi->facet_begin());
771  geomSizes.push_back(cntVertInAgg);
772  }
773  }
774  }
775 
776  }
777 #endif
778 
779  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
780  void VisualizationHelpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::doGraphEdges(std::vector<int>& vertices, std::vector<int>& geomSizes, Teuchos::RCP<GraphBase>& G, Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> & /* fx */, Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> & /* fy */, Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> & /* fz */) {
781  ArrayView<const Scalar> values;
782  ArrayView<const LocalOrdinal> neighbors;
783 
784  std::vector<std::pair<int, int> > vert1; //vertices (node indices)
785 
786  ArrayView<const LocalOrdinal> indices;
787  for(LocalOrdinal locRow = 0; locRow < LocalOrdinal(G->GetNodeNumVertices()); locRow++) {
788  neighbors = G->getNeighborVertices(locRow);
789  //Add those local indices (columns) to the list of connections (which are pairs of the form (localM, localN))
790  for(int gEdge = 0; gEdge < int(neighbors.size()); ++gEdge) {
791  vert1.push_back(std::pair<int, int>(locRow, neighbors[gEdge]));
792  }
793  }
794  for(size_t i = 0; i < vert1.size(); i ++) {
795  if(vert1[i].first > vert1[i].second) {
796  int temp = vert1[i].first;
797  vert1[i].first = vert1[i].second;
798  vert1[i].second = temp;
799  }
800  }
801  std::sort(vert1.begin(), vert1.end());
802  std::vector<std::pair<int, int> >::iterator newEnd = unique(vert1.begin(), vert1.end()); //remove duplicate edges
803  vert1.erase(newEnd, vert1.end());
804  //std::vector<int> points1;
805  vertices.reserve(2 * vert1.size());
806  geomSizes.reserve(vert1.size());
807  for(size_t i = 0; i < vert1.size(); i++) {
808  vertices.push_back(vert1[i].first);
809  vertices.push_back(vert1[i].second);
810  geomSizes.push_back(2);
811  }
812  }
813 
814  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
816  {
817  const double ccw_zero_threshold=1e-8;
818  double val = (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
819  return (val > ccw_zero_threshold) ? 1 : ( (val < -ccw_zero_threshold) ? -1 : 0);
820  }
821 
822  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
824  {
825  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);
826  }
827 
828  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
830  {
831  return v1.x * v2.x + v1.y * v2.y;
832  }
833 
834  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
836  {
837  return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z;
838  }
839 
840  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
842  {
843  myVec3 rel(point.x - inPlane.x, point.y - inPlane.y, point.z - inPlane.z); //position of the point relative to the plane
844  return dotProduct(rel, n) > 1e-12 ? true : false;
845  }
846 
847  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
849  {
850  return sqrt(dotProduct(vec, vec));
851  }
852 
853  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
855  {
856  return sqrt(dotProduct(vec, vec));
857  }
858 
859  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
861  {
862  return mymagnitude(vecSubtract(p1, p2));
863  }
864 
865 
866  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
868  {
869  return mymagnitude(vecSubtract(p1, p2));
870  }
871 
872  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
874  {
875  return myVec2(v1.x - v2.x, v1.y - v2.y);
876  }
877 
878  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
880  {
881  return myVec3(v1.x - v2.x, v1.y - v2.y, v1.z - v2.z);
882  }
883 
884  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
885  myVec2 VisualizationHelpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getNorm(myVec2 v) //"normal" to a 2D vector - just rotate 90 degrees to left
886  {
887  return myVec2(v.y, -v.x);
888  }
889 
890  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
891  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)
892  {
893  return crossProduct(vecSubtract(v2, v1), vecSubtract(v3, v1));
894  }
895 
896  //get minimum distance from 'point' to plane containing v1, v2, v3 (or the triangle with v1, v2, v3 as vertices)
897  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
899  {
900  using namespace std;
901  myVec3 norm = getNorm(v1, v2, v3);
902  //must normalize the normal vector
903  double normScl = mymagnitude(norm);
904  double rv = 0.0;
905  if (normScl > 1e-8) {
906  norm.x /= normScl;
907  norm.y /= normScl;
908  norm.z /= normScl;
909  rv = fabs(dotProduct(norm, vecSubtract(point, v1)));
910  } else {
911  // triangle is degenerated
912  myVec3 test1 = vecSubtract(v3, v1);
913  myVec3 test2 = vecSubtract(v2, v1);
914  bool useTest1 = mymagnitude(test1) > 0.0 ? true : false;
915  bool useTest2 = mymagnitude(test2) > 0.0 ? true : false;
916  if(useTest1 == true) {
917  double normScl1 = mymagnitude(test1);
918  test1.x /= normScl1;
919  test1.y /= normScl1;
920  test1.z /= normScl1;
921  rv = fabs(dotProduct(test1, vecSubtract(point, v1)));
922  } else if (useTest2 == true) {
923  double normScl2 = mymagnitude(test2);
924  test2.x /= normScl2;
925  test2.y /= normScl2;
926  test2.z /= normScl2;
927  rv = fabs(dotProduct(test2, vecSubtract(point, v1)));
928  } else {
929  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError,
930  "VisualizationHelpers::pointDistFromTri: Could not determine the distance of a point to a triangle.");
931  }
932 
933  }
934  return rv;
935  }
936 
937  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
938  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 typename Teuchos::ScalarTraits<Scalar>::coordinateType>& xCoords, const Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType>& yCoords, const Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType>& zCoords) {
939  //*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.
940  //precondition: each triangle is already oriented so that getNorm_(v1, v2, v3) points outward (away from interior of hull)
941  //First find the point furthest from triangle.
942  using namespace std;
943  typedef std::list<int>::iterator Iter;
944  typedef std::list<myTriangle>::iterator TriIter;
945  typedef list<pair<int, int> >::iterator EdgeIter;
946  double maxDist = 0;
947  //Need vector representations of triangle's vertices
948  myVec3 v1(xCoords[tri.v1], yCoords[tri.v1], zCoords[tri.v1]);
949  myVec3 v2(xCoords[tri.v2], yCoords[tri.v2], zCoords[tri.v2]);
950  myVec3 v3(xCoords[tri.v3], yCoords[tri.v3], zCoords[tri.v3]);
951  myVec3 farPointVec; //useful to have both the point's coordinates and it's position in the list
952  Iter farPoint = pointsInFront.begin();
953  for(Iter point = pointsInFront.begin(); point != pointsInFront.end(); point++)
954  {
955  myVec3 pointVec(xCoords[*point], yCoords[*point], zCoords[*point]);
956  double dist = pointDistFromTri(pointVec, v1, v2, v3);
957  // Want to make sure nodeDist is significantly larger than maxDist to prevent roundoff errors from impacting node choice.
958  if(dist > maxDist && dist - maxDist > 1e-12)
959  {
960  dist = maxDist;
961  farPointVec = pointVec;
962  farPoint = point;
963  }
964  }
965  //Find all the triangles that the point is in front of (can be more than 1)
966  //At the same time, remove them from tris, as every one will be replaced later
967  vector<myTriangle> visible; //use a list of iterators so that the underlying object is still in tris
968  for(TriIter it = tris.begin(); it != tris.end();)
969  {
970  myVec3 vec1(xCoords[it->v1], yCoords[it->v1], zCoords[it->v1]);
971  myVec3 vec2(xCoords[it->v2], yCoords[it->v2], zCoords[it->v2]);
972  myVec3 vec3(xCoords[it->v3], yCoords[it->v3], zCoords[it->v3]);
973  myVec3 norm = getNorm(vec1, vec2, vec3);
974  if(isInFront(farPointVec, vec1, norm))
975  {
976  visible.push_back(*it);
977  it = tris.erase(it);
978  }
979  else
980  it++;
981  }
982  //Figure out what triangles need to be destroyed/created
983  //First create a list of edges (as std::pair<int, int>, where the two ints are the node endpoints)
984  list<pair<int, int> > horizon;
985  //For each triangle, add edges to the list iff the edge only appears once in the set of all
986  //Have members of horizon have the lower node # first, and the higher one second
987  for(vector<myTriangle>::iterator it = visible.begin(); it != visible.end(); it++)
988  {
989  pair<int, int> e1(it->v1, it->v2);
990  pair<int, int> e2(it->v2, it->v3);
991  pair<int, int> e3(it->v1, it->v3);
992  //"sort" the pair values
993  if(e1.first > e1.second)
994  {
995  int temp = e1.first;
996  e1.first = e1.second;
997  e1.second = temp;
998  }
999  if(e2.first > e2.second)
1000  {
1001  int temp = e2.first;
1002  e2.first = e2.second;
1003  e2.second = temp;
1004  }
1005  if(e3.first > e3.second)
1006  {
1007  int temp = e3.first;
1008  e3.first = e3.second;
1009  e3.second = temp;
1010  }
1011  horizon.push_back(e1);
1012  horizon.push_back(e2);
1013  horizon.push_back(e3);
1014  }
1015  //sort based on lower node first, then higher node (lexicographical ordering built in to pair)
1016  horizon.sort();
1017  //Remove all edges from horizon, except those that appear exactly once
1018  {
1019  EdgeIter it = horizon.begin();
1020  while(it != horizon.end())
1021  {
1022  int occur = count(horizon.begin(), horizon.end(), *it);
1023  if(occur > 1)
1024  {
1025  pair<int, int> removeVal = *it;
1026  while(removeVal == *it && !(it == horizon.end()))
1027  it = horizon.erase(it);
1028  }
1029  else
1030  it++;
1031  }
1032  }
1033  //Now make a list of new triangles being created, each of which take 2 vertices from an edge and one from farPoint
1034  list<myTriangle> newTris;
1035  for(EdgeIter it = horizon.begin(); it != horizon.end(); it++)
1036  {
1037  myTriangle t(it->first, it->second, *farPoint);
1038  newTris.push_back(t);
1039  }
1040  //Ensure every new triangle is oriented outwards, using the barycenter of the initial tetrahedron
1041  vector<myTriangle> trisToProcess;
1042  vector<list<int> > newFrontPoints;
1043  for(TriIter it = newTris.begin(); it != newTris.end(); it++)
1044  {
1045  myVec3 t1(xCoords[it->v1], yCoords[it->v1], zCoords[it->v1]);
1046  myVec3 t2(xCoords[it->v2], yCoords[it->v2], zCoords[it->v2]);
1047  myVec3 t3(xCoords[it->v3], yCoords[it->v3], zCoords[it->v3]);
1048  if(isInFront(barycenter, t1, getNorm(t1, t2, t3)))
1049  {
1050  //need to swap two vertices to flip orientation of triangle
1051  int temp = it->v1;
1052  myVec3 tempVec = t1;
1053  it->v1 = it->v2;
1054  t1 = t2;
1055  it->v2 = temp;
1056  t2 = tempVec;
1057  }
1058  myVec3 outwardNorm = getNorm(t1, t2, t3); //now definitely points outwards
1059  //Add the triangle to tris
1060  tris.push_back(*it);
1061  trisToProcess.push_back(tris.back());
1062  //Make a list of the points that are in front of nextToProcess, to be passed in for processing
1063  list<int> newInFront;
1064  for(Iter point = pointsInFront.begin(); point != pointsInFront.end();)
1065  {
1066  myVec3 pointVec(xCoords[*point], yCoords[*point], zCoords[*point]);
1067  if(isInFront(pointVec, t1, outwardNorm))
1068  {
1069  newInFront.push_back(*point);
1070  point = pointsInFront.erase(point);
1071  }
1072  else
1073  point++;
1074  }
1075  newFrontPoints.push_back(newInFront);
1076  }
1077  vector<myTriangle> allRemoved; //list of all invalid iterators that were erased by calls to processmyTriangle below
1078  for(int i = 0; i < int(trisToProcess.size()); i++)
1079  {
1080  if(!newFrontPoints[i].empty())
1081  {
1082  //Comparing the 'triangle to process' to the one for this call prevents infinite recursion/stack overflow.
1083  //TODO: Why was it doing that? Rounding error? Make more robust fix. But this does work for the time being.
1084  if(find(allRemoved.begin(), allRemoved.end(), trisToProcess[i]) == allRemoved.end() && !(trisToProcess[i] == tri))
1085  {
1086  vector<myTriangle> removedList = processTriangle(tris, trisToProcess[i], newFrontPoints[i], barycenter, xCoords, yCoords, zCoords);
1087  for(int j = 0; j < int(removedList.size()); j++)
1088  allRemoved.push_back(removedList[j]);
1089  }
1090  }
1091  }
1092  return visible;
1093  }
1094 
1095  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1096  std::vector<int> VisualizationHelpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::giftWrap(std::vector<myVec2>& points, std::vector<int>& nodes, const Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> & xCoords, const Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> & yCoords) {
1097  TEUCHOS_TEST_FOR_EXCEPTION(points.size() < 3, Exceptions::RuntimeError,
1098  "CoarseningVisualization::giftWrap: Gift wrap algorithm input has to have at least 3 points!");
1099 
1100 #if 1 // TAW's version to determine "minimal" node
1101  // determine minimal x and y coordinates
1102  double min_x =points[0].x;
1103  double min_y =points[0].y;
1104  for(std::vector<int>::iterator it = nodes.begin(); it != nodes.end(); it++) {
1105  int i = it - nodes.begin();
1106  if(points[i].x < min_x) min_x = points[i].x;
1107  if(points[i].y < min_y) min_y = points[i].y;
1108  }
1109  // create dummy min coordinates
1110  min_x -= 1.0;
1111  min_y -= 1.0;
1112  myVec2 dummy_min(min_x, min_y);
1113 
1114  // loop over all nodes and determine nodes with minimal distance to (min_x, min_y)
1115  std::vector<int> hull;
1116  myVec2 min = points[0];
1117  double mindist = distance(min,dummy_min);
1118  std::vector<int>::iterator minNode = nodes.begin();
1119  for(std::vector<int>::iterator it = nodes.begin(); it != nodes.end(); it++) {
1120  int i = it - nodes.begin();
1121  if(distance(points[i],dummy_min) < mindist) {
1122  mindist = distance(points[i],dummy_min);
1123  min = points[i];
1124  minNode = it;
1125  }
1126  }
1127  hull.push_back(*minNode);
1128 #else // Brian's code
1129  std::vector<int> hull;
1130  std::vector<int>::iterator minNode = nodes.begin();
1131  myVec2 min = points[0];
1132  for(std::vector<int>::iterator it = nodes.begin(); it != nodes.end(); it++)
1133  {
1134  int i = it - nodes.begin();
1135  if(points[i].x < min.x || (fabs(points[i].x - min.x) < 1e-10 && points[i].y < min.y))
1136  {
1137  min = points[i];
1138  minNode = it;
1139  }
1140  }
1141  hull.push_back(*minNode);
1142 #endif
1143 
1144  bool includeMin = false;
1145  //int debug_it = 0;
1146  while(1)
1147  {
1148  std::vector<int>::iterator leftMost = nodes.begin();
1149  if(!includeMin && leftMost == minNode)
1150  {
1151  leftMost++;
1152  }
1153  std::vector<int>::iterator it = leftMost;
1154  it++;
1155  for(; it != nodes.end(); it++)
1156  {
1157  if(it == minNode && !includeMin) //don't compare to min on very first sweep
1158  continue;
1159  if(*it == hull.back())
1160  continue;
1161  //see if it is in front of line containing nodes thisHull.back() and leftMost
1162  //first get the left normal of leftMost - thisHull.back() (<dy, -dx>)
1163  myVec2 leftMostVec = points[leftMost - nodes.begin()];
1164  myVec2 lastVec(xCoords[hull.back()], yCoords[hull.back()]);
1165  myVec2 testNorm = getNorm(vecSubtract(leftMostVec, lastVec));
1166  //now dot testNorm with *it - leftMost. If dot is positive, leftMost becomes it. If dot is zero, take one further from thisHull.back().
1167  myVec2 itVec(xCoords[*it], yCoords[*it]);
1168  double dotProd = dotProduct(testNorm, vecSubtract(itVec, lastVec));
1169  if(-1e-8 < dotProd && dotProd < 1e-8)
1170  {
1171  //thisHull.back(), it and leftMost are collinear.
1172  //Just sum the differences in x and differences in y for each and compare to get further one, don't need distance formula
1173  myVec2 itDist = vecSubtract(itVec, lastVec);
1174  myVec2 leftMostDist = vecSubtract(leftMostVec, lastVec);
1175  if(fabs(itDist.x) + fabs(itDist.y) > fabs(leftMostDist.x) + fabs(leftMostDist.y)) {
1176  leftMost = it;
1177  }
1178  }
1179  else if(dotProd > 0) {
1180  leftMost = it;
1181 
1182  }
1183  }
1184  //if leftMost is min, then the loop is complete.
1185  if(*leftMost == *minNode)
1186  break;
1187  hull.push_back(*leftMost);
1188  includeMin = true; //have found second point (the one after min) so now include min in the searches
1189  //debug_it ++;
1190  //if(debug_it > 100) exit(0); //break;
1191  }
1192  return hull;
1193  }
1194 
1195  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1196  std::vector<int> VisualizationHelpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::makeUnique(std::vector<int>& vertices) const
1197  {
1198  using namespace std;
1199  vector<int> uniqueNodes = vertices;
1200  sort(uniqueNodes.begin(), uniqueNodes.end());
1201  vector<int>::iterator newUniqueFineEnd = unique(uniqueNodes.begin(), uniqueNodes.end());
1202  uniqueNodes.erase(newUniqueFineEnd, uniqueNodes.end());
1203  //uniqueNodes is now a sorted list of the nodes whose info actually goes in file
1204  //Now replace values in vertices with locations of the old values in uniqueFine
1205  for(int i = 0; i < int(vertices.size()); i++)
1206  {
1207  int lo = 0;
1208  int hi = uniqueNodes.size() - 1;
1209  int mid = 0;
1210  int search = vertices[i];
1211  while(lo <= hi)
1212  {
1213  mid = lo + (hi - lo) / 2;
1214  if(uniqueNodes[mid] == search)
1215  break;
1216  else if(uniqueNodes[mid] > search)
1217  hi = mid - 1;
1218  else
1219  lo = mid + 1;
1220  }
1221  if(uniqueNodes[mid] != search)
1222  throw runtime_error("Issue in makeUnique_() - a point wasn't found in list.");
1223  vertices[i] = mid;
1224  }
1225  return uniqueNodes;
1226  }
1227 
1228  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1229  std::string VisualizationHelpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::replaceAll(std::string result, const std::string& replaceWhat, const std::string& replaceWithWhat) const {
1230  while(1) {
1231  const int pos = result.find(replaceWhat);
1232  if (pos == -1)
1233  break;
1234  result.replace(pos, replaceWhat.size(), replaceWithWhat);
1235  }
1236  return result;
1237  }
1238 
1239  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1240  std::string VisualizationHelpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getFileName(int numProcs, int myRank, int level, const Teuchos::ParameterList & pL) const {
1241  std::string filenameToWrite = getBaseFileName(numProcs, level, pL);
1242  filenameToWrite = this->replaceAll(filenameToWrite, "%PROCID", toString(myRank));
1243  return filenameToWrite;
1244  }
1245 
1246  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1247  std::string VisualizationHelpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getBaseFileName(int numProcs, int level, const Teuchos::ParameterList & pL) const {
1248  std::string filenameToWrite = pL.get<std::string>("visualization: output filename");
1249  int timeStep = pL.get<int>("visualization: output file: time step");
1250  int iter = pL.get<int>("visualization: output file: iter");
1251 
1252  if(filenameToWrite.rfind(".vtu") == std::string::npos)
1253  filenameToWrite.append(".vtu");
1254  if(numProcs > 1 && filenameToWrite.rfind("%PROCID") == std::string::npos) //filename can't be identical between processsors in parallel problem
1255  filenameToWrite.insert(filenameToWrite.rfind(".vtu"), "-proc%PROCID");
1256 
1257  filenameToWrite = this->replaceAll(filenameToWrite, "%LEVELID", toString(level));
1258  filenameToWrite = this->replaceAll(filenameToWrite, "%TIMESTEP", toString(timeStep));
1259  filenameToWrite = this->replaceAll(filenameToWrite, "%ITER", toString(iter));
1260  return filenameToWrite;
1261  }
1262 
1263  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1264  std::string VisualizationHelpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getPVTUFileName(int numProcs, int /* myRank */, int level, const Teuchos::ParameterList & pL) const {
1265  std::string filenameToWrite = getBaseFileName(numProcs, level, pL);
1266  std::string masterStem = filenameToWrite.substr(0, filenameToWrite.rfind(".vtu"));
1267  masterStem = this->replaceAll(masterStem, "%PROCID", "");
1268  std::string pvtuFilename = masterStem + "-master.pvtu";
1269  return pvtuFilename;
1270  }
1271 
1272  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1273  void VisualizationHelpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeFileVTKOpening(std::ofstream & fout, std::vector<int> & uniqueFine, std::vector<int> & geomSizesFine) const {
1274  std::string styleName = "PointCloud"; // TODO adapt this
1275 
1276  std::string indent = " ";
1277  fout << "<!--" << styleName << " Aggregates Visualization-->" << std::endl;
1278  fout << "<VTKFile type=\"UnstructuredGrid\" byte_order=\"LittleEndian\">" << std::endl;
1279  fout << " <UnstructuredGrid>" << std::endl;
1280  fout << " <Piece NumberOfPoints=\"" << uniqueFine.size() << "\" NumberOfCells=\"" << geomSizesFine.size() << "\">" << std::endl;
1281  fout << " <PointData Scalars=\"Node Aggregate Processor\">" << std::endl;
1282  }
1283 
1284  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1285  void VisualizationHelpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeFileVTKNodes(std::ofstream & fout, std::vector<int> & uniqueFine, Teuchos::RCP<const Map> & nodeMap) const {
1286  std::string indent = " ";
1287  fout << " <DataArray type=\"Int32\" Name=\"Node\" format=\"ascii\">" << std::endl;
1288  indent = " ";
1289  bool localIsGlobal = GlobalOrdinal(nodeMap->getGlobalNumElements()) == GlobalOrdinal(nodeMap->getNodeNumElements());
1290  for(size_t i = 0; i < uniqueFine.size(); i++)
1291  {
1292  if(localIsGlobal)
1293  {
1294  fout << uniqueFine[i] << " "; //if all nodes are on this processor, do not map from local to global
1295  }
1296  else
1297  fout << nodeMap->getGlobalElement(uniqueFine[i]) << " ";
1298  if(i % 10 == 9)
1299  fout << std::endl << indent;
1300  }
1301  fout << std::endl;
1302  fout << " </DataArray>" << std::endl;
1303  }
1304 
1305  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1306  void VisualizationHelpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeFileVTKData(std::ofstream & fout, std::vector<int> & uniqueFine, LocalOrdinal myAggOffset, ArrayRCP<LocalOrdinal> & vertex2AggId, int myRank) const {
1307  std::string indent = " ";
1308  fout << " <DataArray type=\"Int32\" Name=\"Aggregate\" format=\"ascii\">" << std::endl;
1309  fout << indent;
1310  for(int i = 0; i < int(uniqueFine.size()); i++)
1311  {
1312  fout << myAggOffset + vertex2AggId[uniqueFine[i]] << " ";
1313  if(i % 10 == 9)
1314  fout << std::endl << indent;
1315  }
1316  fout << std::endl;
1317  fout << " </DataArray>" << std::endl;
1318  fout << " <DataArray type=\"Int32\" Name=\"Processor\" format=\"ascii\">" << std::endl;
1319  fout << indent;
1320  for(int i = 0; i < int(uniqueFine.size()); i++)
1321  {
1322  fout << myRank << " ";
1323  if(i % 20 == 19)
1324  fout << std::endl << indent;
1325  }
1326  fout << std::endl;
1327  fout << " </DataArray>" << std::endl;
1328  fout << " </PointData>" << std::endl;
1329  }
1330 
1331  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1332  void VisualizationHelpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeFileVTKCoordinates(std::ofstream & fout, std::vector<int> & uniqueFine, Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> & fx, Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> & fy, Teuchos::ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> & fz, int dim) const {
1333  std::string indent = " ";
1334  fout << " <Points>" << std::endl;
1335  fout << " <DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"ascii\">" << std::endl;
1336  fout << indent;
1337  for(int i = 0; i < int(uniqueFine.size()); i++)
1338  {
1339  fout << fx[uniqueFine[i]] << " " << fy[uniqueFine[i]] << " ";
1340  if(dim == 2)
1341  fout << "0 ";
1342  else
1343  fout << fz[uniqueFine[i]] << " ";
1344  if(i % 3 == 2)
1345  fout << std::endl << indent;
1346  }
1347  fout << std::endl;
1348  fout << " </DataArray>" << std::endl;
1349  fout << " </Points>" << std::endl;
1350  }
1351 
1352  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1353  void VisualizationHelpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeFileVTKCells(std::ofstream & fout, std::vector<int> & /* uniqueFine */, std::vector<LocalOrdinal> & vertices, std::vector<LocalOrdinal> & geomSize) const {
1354  std::string indent = " ";
1355  fout << " <Cells>" << std::endl;
1356  fout << " <DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">" << std::endl;
1357  fout << indent;
1358  for(int i = 0; i < int(vertices.size()); i++)
1359  {
1360  fout << vertices[i] << " ";
1361  if(i % 10 == 9)
1362  fout << std::endl << indent;
1363  }
1364  fout << std::endl;
1365  fout << " </DataArray>" << std::endl;
1366  fout << " <DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">" << std::endl;
1367  fout << indent;
1368  int accum = 0;
1369  for(int i = 0; i < int(geomSize.size()); i++)
1370  {
1371  accum += geomSize[i];
1372  fout << accum << " ";
1373  if(i % 10 == 9)
1374  fout << std::endl << indent;
1375  }
1376  fout << std::endl;
1377  fout << " </DataArray>" << std::endl;
1378  fout << " <DataArray type=\"Int32\" Name=\"types\" format=\"ascii\">" << std::endl;
1379  fout << indent;
1380  for(int i = 0; i < int(geomSize.size()); i++)
1381  {
1382  switch(geomSize[i])
1383  {
1384  case 1:
1385  fout << "1 "; //Point
1386  break;
1387  case 2:
1388  fout << "3 "; //Line
1389  break;
1390  case 3:
1391  fout << "5 "; //Triangle
1392  break;
1393  default:
1394  fout << "7 "; //Polygon
1395  }
1396  if(i % 30 == 29)
1397  fout << std::endl << indent;
1398  }
1399  fout << std::endl;
1400  fout << " </DataArray>" << std::endl;
1401  fout << " </Cells>" << std::endl;
1402  }
1403 
1404  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1406  fout << " </Piece>" << std::endl;
1407  fout << " </UnstructuredGrid>" << std::endl;
1408  fout << "</VTKFile>" << std::endl;
1409  }
1410 
1411 
1412  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1413  void VisualizationHelpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writePVTU(std::ofstream& pvtu, std::string baseFname, int numProcs, bool bFineEdges, bool /* bCoarseEdges */) const {
1414  //If using vtk, filenameToWrite now contains final, correct ***.vtu filename (for the current proc)
1415  //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
1416  //pvtu file.
1417  pvtu << "<VTKFile type=\"PUnstructuredGrid\" byte_order=\"LittleEndian\">" << std::endl;
1418  pvtu << " <PUnstructuredGrid GhostLevel=\"0\">" << std::endl;
1419  pvtu << " <PPointData Scalars=\"Node Aggregate Processor\">" << std::endl;
1420  pvtu << " <PDataArray type=\"Int32\" Name=\"Node\"/>" << std::endl;
1421  pvtu << " <PDataArray type=\"Int32\" Name=\"Aggregate\"/>" << std::endl;
1422  pvtu << " <PDataArray type=\"Int32\" Name=\"Processor\"/>" << std::endl;
1423  pvtu << " </PPointData>" << std::endl;
1424  pvtu << " <PPoints>" << std::endl;
1425  pvtu << " <PDataArray type=\"Float64\" NumberOfComponents=\"3\"/>" << std::endl;
1426  pvtu << " </PPoints>" << std::endl;
1427  for(int i = 0; i < numProcs; i++) {
1428  //specify the piece for each proc (the replaceAll expression matches what the filenames will be on other procs)
1429  pvtu << " <Piece Source=\"" << replaceAll(baseFname, "%PROCID", toString(i)) << "\"/>" << std::endl;
1430  }
1431  //reference files with graph pieces, if applicable
1432  if(bFineEdges)
1433  {
1434  for(int i = 0; i < numProcs; i++)
1435  {
1436  std::string fn = replaceAll(baseFname, "%PROCID", toString(i));
1437  pvtu << " <Piece Source=\"" << fn.insert(fn.rfind(".vtu"), "-finegraph") << "\"/>" << std::endl;
1438  }
1439  }
1440  /*if(doCoarseGraphEdges_)
1441  {
1442  for(int i = 0; i < numProcs; i++)
1443  {
1444  std::string fn = replaceAll(baseFname, "%PROCID", toString(i));
1445  pvtu << " <Piece Source=\"" << fn.insert(fn.rfind(".vtu"), "-coarsegraph") << "\"/>" << std::endl;
1446  }
1447  }*/
1448  pvtu << " </PUnstructuredGrid>" << std::endl;
1449  pvtu << "</VTKFile>" << std::endl;
1450  pvtu.close();
1451  }
1452 
1453  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1455  try {
1456  std::ofstream color("random-colormap.xml");
1457  color << "<ColorMap name=\"MueLu-Random\" space=\"RGB\">" << std::endl;
1458  //Give -1, -2, -3 distinctive colors (so that the style functions can have constrasted geometry types)
1459  //Do red, orange, yellow to constrast with cool color spectrum for other types
1460  color << " <Point x=\"" << -1 << "\" o=\"1\" r=\"1\" g=\"0\" b=\"0\"/>" << std::endl;
1461  color << " <Point x=\"" << -2 << "\" o=\"1\" r=\"1\" g=\"0.6\" b=\"0\"/>" << std::endl;
1462  color << " <Point x=\"" << -3 << "\" o=\"1\" r=\"1\" g=\"1\" b=\"0\"/>" << std::endl;
1463  srand(time(NULL));
1464  for(int i = 0; i < 5000; i += 4) {
1465  color << " <Point x=\"" << i << "\" o=\"1\" r=\"" << (rand() % 50) / 256.0 << "\" g=\"" << (rand() % 256) / 256.0 << "\" b=\"" << (rand() % 256) / 256.0 << "\"/>" << std::endl;
1466  }
1467  color << "</ColorMap>" << std::endl;
1468  color.close();
1469  }
1470  catch(std::exception& e) {
1471  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError,
1472  "VisualizationHelpers::buildColormap: Error while building colormap file: " << e.what());
1473  }
1474  }
1475 
1476 } // namespace MueLu
1477 
1478 #endif /* MUELU_VISUALIZATIONHELPERS_DEF_HPP_ */
void writePVTU(std::ofstream &pvtu, std::string baseFname, int numProcs, bool bFineEdges=false, bool bCoarseEdges=false) const
static void doPointCloud(std::vector< int > &vertices, std::vector< int > &geomSizes, LO numLocalAggs, LO numFineNodes)
void writeFileVTKOpening(std::ofstream &fout, std::vector< int > &uniqueFine, std::vector< int > &geomSizesFine) const
MueLu::DefaultLocalOrdinal LocalOrdinal
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
RCP< ParameterList > GetValidParameterList() const
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 typename Teuchos::ScalarTraits< Scalar >::coordinateType > &xCoords, const Teuchos::ArrayRCP< const typename Teuchos::ScalarTraits< Scalar >::coordinateType > &yCoords)
static void doCGALConvexHulls3D(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 typename Teuchos::ScalarTraits< Scalar >::coordinateType > &xCoords, const Teuchos::ArrayRCP< const typename Teuchos::ScalarTraits< Scalar >::coordinateType > &yCoords, const Teuchos::ArrayRCP< const typename Teuchos::ScalarTraits< Scalar >::coordinateType > &zCoords)
void writeFileVTKCoordinates(std::ofstream &fout, std::vector< int > &uniqueFine, Teuchos::ArrayRCP< const typename Teuchos::ScalarTraits< Scalar >::coordinateType > &fx, Teuchos::ArrayRCP< const typename Teuchos::ScalarTraits< Scalar >::coordinateType > &fy, Teuchos::ArrayRCP< const typename Teuchos::ScalarTraits< Scalar >::coordinateType > &fz, int dim) const
static void doCGALConvexHulls2D(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 typename Teuchos::ScalarTraits< Scalar >::coordinateType > &xCoords, const Teuchos::ArrayRCP< const typename Teuchos::ScalarTraits< Scalar >::coordinateType > &yCoords)
std::string getPVTUFileName(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
static int ccw(const myVec2 &a, const myVec2 &b, const myVec2 &c)
static double pointDistFromTri(myVec3 point, myVec3 v1, myVec3 v2, myVec3 v3)
static double dotProduct(myVec2 v1, myVec2 v2)
std::vector< int > makeUnique(std::vector< int > &vertices) const
replaces node indices in vertices with compressed unique indices, and returns list of unique points ...
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 typename Teuchos::ScalarTraits< Scalar >::coordinateType > &xCoords, const Teuchos::ArrayRCP< const typename Teuchos::ScalarTraits< Scalar >::coordinateType > &yCoords, const Teuchos::ArrayRCP< const typename Teuchos::ScalarTraits< Scalar >::coordinateType > &zCoords)
MueLu::DefaultGlobalOrdinal GlobalOrdinal
std::string replaceAll(std::string result, const std::string &replaceWhat, const std::string &replaceWithWhat) const
static void doGraphEdges(std::vector< int > &vertices, std::vector< int > &geomSizes, Teuchos::RCP< GraphBase > &G, Teuchos::ArrayRCP< const typename Teuchos::ScalarTraits< Scalar >::coordinateType > &fx, Teuchos::ArrayRCP< const typename Teuchos::ScalarTraits< Scalar >::coordinateType > &fy, Teuchos::ArrayRCP< const typename Teuchos::ScalarTraits< Scalar >::coordinateType > &fz)
static void doJacks(std::vector< int > &vertices, std::vector< int > &geomSizes, LO numLocalAggs, LO numFineNodes, const std::vector< bool > &isRoot, const ArrayRCP< LO > &vertex2AggId)
void writeFileVTKCells(std::ofstream &fout, std::vector< int > &uniqueFine, std::vector< LocalOrdinal > &vertices, std::vector< LocalOrdinal > &geomSize) const
void writeFileVTKData(std::ofstream &fout, std::vector< int > &uniqueFine, LocalOrdinal myAggOffset, ArrayRCP< LocalOrdinal > &vertex2AggId, int myRank) const
static myVec3 crossProduct(myVec3 v1, myVec3 v2)
static myVec2 vecSubtract(myVec2 v1, myVec2 v2)
std::string getFileName(int numProcs, int myRank, int level, const Teuchos::ParameterList &pL) const
void writeFileVTKClosing(std::ofstream &fout) const
static double distance(myVec2 p1, myVec2 p2)
static bool isInFront(myVec3 point, myVec3 inPlane, myVec3 n)
Exception throws to report errors in the internal logical of the program.
std::string getBaseFileName(int numProcs, int level, const Teuchos::ParameterList &pL) const
static std::vector< myTriangle > processTriangle(std::list< myTriangle > &tris, myTriangle tri, std::list< int > &pointsInFront, myVec3 &barycenter, const Teuchos::ArrayRCP< const typename Teuchos::ScalarTraits< Scalar >::coordinateType > &xCoords, const Teuchos::ArrayRCP< const typename Teuchos::ScalarTraits< Scalar >::coordinateType > &yCoords, const Teuchos::ArrayRCP< const typename Teuchos::ScalarTraits< Scalar >::coordinateType > &zCoords)
void replaceAll(std::string &str, const std::string &from, const std::string &to)
static std::vector< int > giftWrap(std::vector< myVec2 > &points, std::vector< int > &nodes, const Teuchos::ArrayRCP< const typename Teuchos::ScalarTraits< Scalar >::coordinateType > &xCoords, const Teuchos::ArrayRCP< const typename Teuchos::ScalarTraits< Scalar >::coordinateType > &yCoords)