MueLu Version of the Day
Loading...
Searching...
No Matches
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
71namespace 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) {
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);
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);
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 });
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>
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->getLocalNumElements());
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_ */
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Exception throws to report errors in the internal logical of the program.
void writeFileVTKData(std::ofstream &fout, std::vector< int > &uniqueFine, LocalOrdinal myAggOffset, ArrayRCP< LocalOrdinal > &vertex2AggId, int myRank) const
static int ccw(const myVec2 &a, const myVec2 &b, const myVec2 &c)
static double pointDistFromTri(myVec3 point, myVec3 v1, myVec3 v2, myVec3 v3)
void writeFileVTKClosing(std::ofstream &fout) const
static myVec3 crossProduct(myVec3 v1, myVec3 v2)
std::string replaceAll(std::string result, const std::string &replaceWhat, const std::string &replaceWithWhat) const
static myVec2 vecSubtract(myVec2 v1, myVec2 v2)
std::string getBaseFileName(int numProcs, int level, const Teuchos::ParameterList &pL) const
std::string getFileName(int numProcs, int myRank, int level, const Teuchos::ParameterList &pL) const
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::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 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)
static double dotProduct(myVec2 v1, myVec2 v2)
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 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)
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)
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 bool isInFront(myVec3 point, myVec3 inPlane, myVec3 n)
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 void doPointCloud(std::vector< int > &vertices, std::vector< int > &geomSizes, LO numLocalAggs, LO numFineNodes)
RCP< ParameterList > GetValidParameterList() const
void writePVTU(std::ofstream &pvtu, std::string baseFname, int numProcs, bool bFineEdges=false, bool bCoarseEdges=false) const
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)
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 writeFileVTKNodes(std::ofstream &fout, std::vector< int > &uniqueFine, Teuchos::RCP< const Map > &nodeMap) 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)
Namespace for MueLu classes and methods.
void replaceAll(std::string &str, const std::string &from, const std::string &to)
std::string toString(const T &what)
Little helper function to convert non-string types to strings.