47 #ifndef MUELU_VISUALIZATIONHELPERS_DEF_HPP_
48 #define MUELU_VISUALIZATIONHELPERS_DEF_HPP_
51 #include <Xpetra_CrsMatrixWrap.hpp>
53 #include <Xpetra_MultiVectorFactory.hpp>
64 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
68 validParamList->
set<std::string>(
"visualization: output filename",
"viz%LEVELID",
"filename for VTK-style visualization output");
69 validParamList->
set<
int>(
"visualization: output file: time step", 0,
"time step variable for output file name");
70 validParamList->
set<
int>(
"visualization: output file: iter", 0,
"nonlinear iteration variable for output file name");
71 validParamList->
set<std::string>(
"visualization: style",
"Point Cloud",
"style of aggregate visualization for VTK output. Can be 'Point Cloud', 'Jacks', 'Convex Hulls'");
72 validParamList->
set<
bool>(
"visualization: build colormap",
false,
"Whether to build a random color map in a separate xml file.");
73 validParamList->
set<
bool>(
"visualization: fine graph edges",
false,
"Whether to draw all fine node connections along with the aggregates.");
75 return validParamList;
78 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
80 vertices.reserve(numFineNodes);
81 geomSizes.reserve(numFineNodes);
83 vertices.push_back(i);
84 geomSizes.push_back(1);
88 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
92 vertices.reserve(vertices.size() + 3 * (numFineNodes - numLocalAggs));
93 geomSizes.reserve(vertices.size() + 2 * (numFineNodes - numLocalAggs));
95 for (
int i = 0; i < numLocalAggs; i++)
99 int numInAggFound = 0;
100 for (
int j = 0; j < numFineNodes; j++) {
106 if (vertex2AggId[root] == vertex2AggId[j]) {
107 vertices.push_back(root);
108 vertices.push_back(j);
109 geomSizes.push_back(2);
111 vertices.push_back(j);
112 geomSizes.push_back(1);
122 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
123 void VisualizationHelpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::doConvexHulls2D(std::vector<int>& vertices, std::vector<int>& geomSizes,
LO numLocalAggs,
LO numFineNodes,
const std::vector<bool>& ,
const ArrayRCP<LO>& vertex2AggId,
const Teuchos::ArrayRCP<
const typename Teuchos::ScalarTraits<Scalar>::coordinateType>& xCoords,
const Teuchos::ArrayRCP<
const typename Teuchos::ScalarTraits<Scalar>::coordinateType>& yCoords) {
126 for (
int agg = 0; agg < numLocalAggs; agg++) {
127 std::vector<int> aggNodes;
128 for (
int i = 0; i < numFineNodes; i++) {
129 if (vertex2AggId[i] == agg)
130 aggNodes.push_back(i);
135 "CoarseningVisualization::doConvexHulls2D: aggregate contains zero nodes!");
136 if (aggNodes.size() == 1) {
137 vertices.push_back(aggNodes.front());
138 geomSizes.push_back(1);
141 if (aggNodes.size() == 2) {
142 vertices.push_back(aggNodes.front());
143 vertices.push_back(aggNodes.back());
144 geomSizes.push_back(2);
147 int N = (int)aggNodes.size();
148 using MyPair = std::pair<myVec2, int>;
149 std::vector<MyPair> pointsAndIndex(N);
150 for (
int i = 0; i < N; i++) {
151 pointsAndIndex[i] = std::make_pair(
myVec2(xCoords[aggNodes[i]], yCoords[aggNodes[i]]), i);
155 std::sort(pointsAndIndex.begin(), pointsAndIndex.end(), [](
const MyPair& a,
const MyPair& b) {
156 return a.first.x < b.first.x || (a.first.x == b.first.x && a.first.y < b.first.y);
160 bool colinear =
true;
161 for (
int i = 0; i < N; i++) {
162 if (ccw(pointsAndIndex[i].first, pointsAndIndex[(i + 1) % N].first, pointsAndIndex[(i + 2) % N].first) != 0) {
169 vertices.push_back(aggNodes[pointsAndIndex.front().second]);
170 vertices.push_back(aggNodes[pointsAndIndex.back().second]);
171 geomSizes.push_back(2);
173 std::vector<int> hull(2 * N);
177 for (
int i = 0; i < N; i++) {
178 while (count >= 2 && ccw(pointsAndIndex[hull[count - 2]].first, pointsAndIndex[hull[count - 1]].first, pointsAndIndex[i].first) <= 0) {
186 for (
int i = N - 1; i > 0; i--) {
187 while (count >= t && ccw(pointsAndIndex[hull[count - 2]].first, pointsAndIndex[hull[count - 1]].first, pointsAndIndex[i - 1].first) <= 0) {
190 hull[count++] = i - 1;
194 hull.resize(count - 1);
197 for (
int i = 0; i < (int)hull.size(); i++) {
198 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");
202 for (
int i = 0; i < (int)hull.size(); i++) {
203 hull[i] = aggNodes[pointsAndIndex[hull[i]].second];
206 vertices.reserve(vertices.size() + hull.size());
207 for (
size_t i = 0; i < hull.size(); i++) {
208 vertices.push_back(hull[i]);
210 geomSizes.push_back(hull.size());
216 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
217 void VisualizationHelpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::doConvexHulls3D(std::vector<int>& vertices, std::vector<int>& geomSizes,
LO numLocalAggs,
LO numFineNodes,
const std::vector<bool>& ,
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) {
222 typedef std::list<int>::iterator Iter;
223 for (
int agg = 0; agg < numLocalAggs; agg++) {
224 std::list<int> aggNodes;
225 for (
int i = 0; i < numFineNodes; i++) {
226 if (vertex2AggId[i] == agg)
227 aggNodes.push_back(i);
231 "CoarseningVisualization::doConvexHulls3D: aggregate contains zero nodes!");
232 if (aggNodes.size() == 1) {
233 vertices.push_back(aggNodes.front());
234 geomSizes.push_back(1);
236 }
else if (aggNodes.size() == 2) {
237 vertices.push_back(aggNodes.front());
238 vertices.push_back(aggNodes.back());
239 geomSizes.push_back(2);
243 bool areCollinear =
true;
245 Iter it = aggNodes.begin();
246 myVec3 firstVec(xCoords[*it], yCoords[*it], zCoords[*it]);
250 myVec3 secondVec(xCoords[*it], yCoords[*it], zCoords[*it]);
251 comp = vecSubtract(secondVec, firstVec);
254 for (; it != aggNodes.end(); it++) {
255 myVec3 thisVec(xCoords[*it], yCoords[*it], zCoords[*it]);
256 myVec3 cross = crossProduct(vecSubtract(thisVec, firstVec), comp);
257 if (mymagnitude(cross) > 1e-10) {
258 areCollinear =
false;
266 Iter min = aggNodes.begin();
267 Iter max = aggNodes.begin();
268 Iter it = ++aggNodes.begin();
269 for (; it != aggNodes.end(); it++) {
270 if (xCoords[*it] < xCoords[*min])
272 else if (xCoords[*it] == xCoords[*min]) {
273 if (yCoords[*it] < yCoords[*min])
275 else if (yCoords[*it] == yCoords[*min]) {
276 if (zCoords[*it] < zCoords[*min]) min = it;
279 if (xCoords[*it] > xCoords[*max])
281 else if (xCoords[*it] == xCoords[*max]) {
282 if (yCoords[*it] > yCoords[*max])
284 else if (yCoords[*it] == yCoords[*max]) {
285 if (zCoords[*it] > zCoords[*max])
290 vertices.push_back(*min);
291 vertices.push_back(*max);
292 geomSizes.push_back(2);
295 bool areCoplanar =
true;
298 Iter vert = aggNodes.begin();
299 myVec3 v1(xCoords[*vert], yCoords[*vert], zCoords[*vert]);
301 myVec3 v2(xCoords[*vert], yCoords[*vert], zCoords[*vert]);
303 myVec3 v3(xCoords[*vert], yCoords[*vert], zCoords[*vert]);
306 while (mymagnitude(crossProduct(vecSubtract(v1, v2), vecSubtract(v1, v3))) < 1e-10) {
308 v3 =
myVec3(xCoords[*vert], yCoords[*vert], zCoords[*vert]);
311 for (; vert != aggNodes.end(); vert++) {
312 myVec3 pt(xCoords[*vert], yCoords[*vert], zCoords[*vert]);
313 if (fabs(pointDistFromTri(pt, v1, v2, v3)) > 1e-12) {
320 myVec3 planeNorm = getNorm(v1, v2, v3);
321 planeNorm.
x = fabs(planeNorm.
x);
322 planeNorm.
y = fabs(planeNorm.
y);
323 planeNorm.
z = fabs(planeNorm.
z);
324 std::vector<myVec2> points;
325 std::vector<int> nodes;
326 if (planeNorm.
x >= planeNorm.
y && planeNorm.
x >= planeNorm.
z) {
328 for (Iter it = aggNodes.begin(); it != aggNodes.end(); it++) {
329 nodes.push_back(*it);
330 points.push_back(
myVec2(yCoords[*it], zCoords[*it]));
333 if (planeNorm.
y >= planeNorm.
x && planeNorm.
y >= planeNorm.
z) {
335 for (Iter it = aggNodes.begin(); it != aggNodes.end(); it++) {
336 nodes.push_back(*it);
337 points.push_back(
myVec2(xCoords[*it], zCoords[*it]));
340 if (planeNorm.
z >= planeNorm.
x && planeNorm.
z >= planeNorm.
y) {
341 for (Iter it = aggNodes.begin(); it != aggNodes.end(); it++) {
342 nodes.push_back(*it);
343 points.push_back(
myVec2(xCoords[*it], yCoords[*it]));
346 std::vector<int> convhull2d = giftWrap(points, nodes, xCoords, yCoords);
347 geomSizes.push_back(convhull2d.size());
348 vertices.reserve(vertices.size() + convhull2d.size());
349 for (
size_t i = 0; i < convhull2d.size(); i++)
350 vertices.push_back(convhull2d[i]);
354 Iter exIt = aggNodes.begin();
355 int extremeSix[] = {*exIt, *exIt, *exIt, *exIt, *exIt, *exIt};
357 for (; exIt != aggNodes.end(); exIt++) {
359 if (xCoords[*it] < xCoords[extremeSix[0]] ||
360 (xCoords[*it] == xCoords[extremeSix[0]] && yCoords[*it] < yCoords[extremeSix[0]]) ||
361 (xCoords[*it] == xCoords[extremeSix[0]] && yCoords[*it] == yCoords[extremeSix[0]] && zCoords[*it] < zCoords[extremeSix[0]]))
363 if (xCoords[*it] > xCoords[extremeSix[1]] ||
364 (xCoords[*it] == xCoords[extremeSix[1]] && yCoords[*it] > yCoords[extremeSix[1]]) ||
365 (xCoords[*it] == xCoords[extremeSix[1]] && yCoords[*it] == yCoords[extremeSix[1]] && zCoords[*it] > zCoords[extremeSix[1]]))
367 if (yCoords[*it] < yCoords[extremeSix[2]] ||
368 (yCoords[*it] == yCoords[extremeSix[2]] && zCoords[*it] < zCoords[extremeSix[2]]) ||
369 (yCoords[*it] == yCoords[extremeSix[2]] && zCoords[*it] == zCoords[extremeSix[2]] && xCoords[*it] < xCoords[extremeSix[2]]))
371 if (yCoords[*it] > yCoords[extremeSix[3]] ||
372 (yCoords[*it] == yCoords[extremeSix[3]] && zCoords[*it] > zCoords[extremeSix[3]]) ||
373 (yCoords[*it] == yCoords[extremeSix[3]] && zCoords[*it] == zCoords[extremeSix[3]] && xCoords[*it] > xCoords[extremeSix[3]]))
375 if (zCoords[*it] < zCoords[extremeSix[4]] ||
376 (zCoords[*it] == zCoords[extremeSix[4]] && xCoords[*it] < xCoords[extremeSix[4]]) ||
377 (zCoords[*it] == zCoords[extremeSix[4]] && xCoords[*it] == xCoords[extremeSix[4]] && yCoords[*it] < yCoords[extremeSix[4]]))
379 if (zCoords[*it] > zCoords[extremeSix[5]] ||
380 (zCoords[*it] == zCoords[extremeSix[5]] && xCoords[*it] > xCoords[extremeSix[5]]) ||
381 (zCoords[*it] == zCoords[extremeSix[5]] && xCoords[*it] == xCoords[extremeSix[5]] && yCoords[*it] > zCoords[extremeSix[5]]))
385 for (
int i = 0; i < 6; i++) {
386 myVec3 thisExtremeVec(xCoords[extremeSix[i]], yCoords[extremeSix[i]], zCoords[extremeSix[i]]);
387 extremeVectors[i] = thisExtremeVec;
392 for (
int i = 0; i < 5; i++) {
393 for (
int j = i + 1; j < 6; j++) {
394 double thisDist = distance(extremeVectors[i], extremeVectors[j]);
396 if (thisDist > maxDist && thisDist - maxDist > 1e-12) {
403 std::list<myTriangle> hullBuilding;
405 aggNodes.remove(extremeSix[base1]);
406 aggNodes.remove(extremeSix[base2]);
409 tri.
v1 = extremeSix[base1];
410 tri.
v2 = extremeSix[base2];
414 myVec3 b1 = extremeVectors[base1];
415 myVec3 b2 = extremeVectors[base2];
417 for (Iter node = aggNodes.begin(); node != aggNodes.end(); node++) {
418 myVec3 nodePos(xCoords[*node], yCoords[*node], zCoords[*node]);
419 double dist = mymagnitude(crossProduct(vecSubtract(nodePos, b1), vecSubtract(nodePos, b2))) / mymagnitude(vecSubtract(b2, b1));
421 if (dist > maxDist && dist - maxDist > 1e-12) {
428 hullBuilding.push_back(tri);
429 myVec3 b3(xCoords[*thirdNode], yCoords[*thirdNode], zCoords[*thirdNode]);
430 aggNodes.erase(thirdNode);
433 int fourthVertex = -1;
434 for (Iter node = aggNodes.begin(); node != aggNodes.end(); node++) {
435 myVec3 thisNode(xCoords[*node], yCoords[*node], zCoords[*node]);
436 double nodeDist = pointDistFromTri(thisNode, b1, b2, b3);
438 if (nodeDist > maxDist && nodeDist - maxDist > 1e-12) {
440 fourthVertex = *node;
443 aggNodes.remove(fourthVertex);
444 myVec3 b4(xCoords[fourthVertex], yCoords[fourthVertex], zCoords[fourthVertex]);
447 tri = hullBuilding.front();
448 tri.
v1 = fourthVertex;
449 hullBuilding.push_back(tri);
450 tri = hullBuilding.front();
451 tri.
v2 = fourthVertex;
452 hullBuilding.push_back(tri);
453 tri = hullBuilding.front();
454 tri.
v3 = fourthVertex;
455 hullBuilding.push_back(tri);
457 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);
458 for (std::list<myTriangle>::iterator tetTri = hullBuilding.begin(); tetTri != hullBuilding.end(); tetTri++) {
459 myVec3 triVert1(xCoords[tetTri->v1], yCoords[tetTri->v1], zCoords[tetTri->v1]);
460 myVec3 triVert2(xCoords[tetTri->v2], yCoords[tetTri->v2], zCoords[tetTri->v2]);
461 myVec3 triVert3(xCoords[tetTri->v3], yCoords[tetTri->v3], zCoords[tetTri->v3]);
462 myVec3 trinorm = getNorm(triVert1, triVert2, triVert3);
463 myVec3 ptInPlane = tetTri == hullBuilding.begin() ? b1 : b4;
464 if (isInFront(barycenter, ptInPlane, trinorm)) {
467 int temp = tetTri->v1;
468 tetTri->v1 = tetTri->v2;
480 for (std::list<myTriangle>::iterator tetTri = hullBuilding.begin(); tetTri != hullBuilding.end(); tetTri++) {
481 myVec3 triVert1(xCoords[tetTri->v1], yCoords[tetTri->v1], zCoords[tetTri->v1]);
482 myVec3 triVert2(xCoords[tetTri->v2], yCoords[tetTri->v2], zCoords[tetTri->v2]);
483 myVec3 triVert3(xCoords[tetTri->v3], yCoords[tetTri->v3], zCoords[tetTri->v3]);
484 trinorms[index] = getNorm(triVert1, triVert2, triVert3);
487 std::list<int> startPoints1;
488 std::list<int> startPoints2;
489 std::list<int> startPoints3;
490 std::list<int> startPoints4;
493 Iter point = aggNodes.begin();
494 while (!aggNodes.empty())
496 myVec3 pointVec(xCoords[*point], yCoords[*point], zCoords[*point]);
499 if (isInFront(pointVec, b1, trinorms[0])) {
500 startPoints1.push_back(*point);
501 point = aggNodes.erase(point);
502 }
else if (isInFront(pointVec, b4, trinorms[1])) {
503 startPoints2.push_back(*point);
504 point = aggNodes.erase(point);
505 }
else if (isInFront(pointVec, b4, trinorms[2])) {
506 startPoints3.push_back(*point);
507 point = aggNodes.erase(point);
508 }
else if (isInFront(pointVec, b4, trinorms[3])) {
509 startPoints4.push_back(*point);
510 point = aggNodes.erase(point);
512 point = aggNodes.erase(point);
517 typedef std::list<myTriangle>::iterator TriIter;
518 TriIter firstTri = hullBuilding.begin();
527 if (!startPoints1.empty())
528 processTriangle(hullBuilding, start1, startPoints1, barycenter, xCoords, yCoords, zCoords);
529 if (!startPoints2.empty())
530 processTriangle(hullBuilding, start2, startPoints2, barycenter, xCoords, yCoords, zCoords);
531 if (!startPoints3.empty())
532 processTriangle(hullBuilding, start3, startPoints3, barycenter, xCoords, yCoords, zCoords);
533 if (!startPoints4.empty())
534 processTriangle(hullBuilding, start4, startPoints4, barycenter, xCoords, yCoords, zCoords);
537 vertices.reserve(vertices.size() + 3 * hullBuilding.size());
538 for (TriIter hullTri = hullBuilding.begin(); hullTri != hullBuilding.end(); hullTri++) {
539 vertices.push_back(hullTri->v1);
540 vertices.push_back(hullTri->v2);
541 vertices.push_back(hullTri->v3);
542 geomSizes.push_back(3);
547 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
551 std::vector<std::pair<int, int> > vert1;
557 for (
int gEdge = 0; gEdge < int(neighbors.length); ++gEdge) {
558 vert1.push_back(std::pair<int, int>(locRow, neighbors(gEdge)));
561 for (
size_t i = 0; i < vert1.size(); i++) {
562 if (vert1[i].first > vert1[i].second) {
563 int temp = vert1[i].first;
564 vert1[i].first = vert1[i].second;
565 vert1[i].second = temp;
568 std::sort(vert1.begin(), vert1.end());
569 std::vector<std::pair<int, int> >::iterator newEnd = unique(vert1.begin(), vert1.end());
570 vert1.erase(newEnd, vert1.end());
572 vertices.reserve(2 * vert1.size());
573 geomSizes.reserve(vert1.size());
574 for (
size_t i = 0; i < vert1.size(); i++) {
575 vertices.push_back(vert1[i].first);
576 vertices.push_back(vert1[i].second);
577 geomSizes.push_back(2);
581 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
583 const double ccw_zero_threshold = 1e-8;
584 double val = (b.
x - a.
x) * (c.
y - a.
y) - (b.
y - a.
y) * (c.
x - a.
x);
585 return (val > ccw_zero_threshold) ? 1 : ((val < -ccw_zero_threshold) ? -1 : 0);
588 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
590 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);
593 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
595 return v1.
x * v2.
x + v1.
y * v2.
y;
598 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
600 return v1.
x * v2.
x + v1.
y * v2.
y + v1.
z * v2.
z;
603 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
605 myVec3 rel(point.
x - inPlane.
x, point.
y - inPlane.
y, point.
z - inPlane.
z);
606 return dotProduct(rel, n) > 1e-12 ?
true :
false;
609 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
611 return sqrt(dotProduct(vec, vec));
614 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
616 return sqrt(dotProduct(vec, vec));
619 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
621 return mymagnitude(vecSubtract(p1, p2));
624 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
626 return mymagnitude(vecSubtract(p1, p2));
629 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
634 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
639 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
645 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
648 return crossProduct(vecSubtract(v2, v1), vecSubtract(v3, v1));
652 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
655 myVec3 norm = getNorm(v1, v2, v3);
657 double normScl = mymagnitude(norm);
659 if (normScl > 1e-8) {
663 rv = fabs(dotProduct(norm, vecSubtract(point, v1)));
666 myVec3 test1 = vecSubtract(v3, v1);
667 myVec3 test2 = vecSubtract(v2, v1);
668 bool useTest1 = mymagnitude(test1) > 0.0 ?
true :
false;
669 bool useTest2 = mymagnitude(test2) > 0.0 ?
true :
false;
670 if (useTest1 ==
true) {
671 double normScl1 = mymagnitude(test1);
675 rv = fabs(dotProduct(test1, vecSubtract(point, v1)));
676 }
else if (useTest2 ==
true) {
677 double normScl2 = mymagnitude(test2);
681 rv = fabs(dotProduct(test2, vecSubtract(point, v1)));
684 "VisualizationHelpers::pointDistFromTri: Could not determine the distance of a point to a triangle.");
690 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
696 typedef std::list<int>::iterator Iter;
697 typedef std::list<myTriangle>::iterator TriIter;
698 typedef list<pair<int, int> >::iterator EdgeIter;
701 myVec3 v1(xCoords[tri.
v1], yCoords[tri.
v1], zCoords[tri.
v1]);
702 myVec3 v2(xCoords[tri.
v2], yCoords[tri.
v2], zCoords[tri.
v2]);
703 myVec3 v3(xCoords[tri.
v3], yCoords[tri.
v3], zCoords[tri.
v3]);
705 Iter farPoint = pointsInFront.begin();
706 for (Iter point = pointsInFront.begin(); point != pointsInFront.end(); point++) {
707 myVec3 pointVec(xCoords[*point], yCoords[*point], zCoords[*point]);
708 double dist = pointDistFromTri(pointVec, v1, v2, v3);
710 if (dist > maxDist && dist - maxDist > 1e-12) {
712 farPointVec = pointVec;
718 vector<myTriangle> visible;
719 for (TriIter it = tris.begin(); it != tris.end();) {
720 myVec3 vec1(xCoords[it->v1], yCoords[it->v1], zCoords[it->v1]);
721 myVec3 vec2(xCoords[it->v2], yCoords[it->v2], zCoords[it->v2]);
722 myVec3 vec3(xCoords[it->v3], yCoords[it->v3], zCoords[it->v3]);
723 myVec3 norm = getNorm(vec1, vec2, vec3);
724 if (isInFront(farPointVec, vec1, norm)) {
725 visible.push_back(*it);
732 list<pair<int, int> > horizon;
735 for (vector<myTriangle>::iterator it = visible.begin(); it != visible.end(); it++) {
736 pair<int, int> e1(it->v1, it->v2);
737 pair<int, int> e2(it->v2, it->v3);
738 pair<int, int> e3(it->v1, it->v3);
740 if (e1.first > e1.second) {
742 e1.first = e1.second;
745 if (e2.first > e2.second) {
747 e2.first = e2.second;
750 if (e3.first > e3.second) {
752 e3.first = e3.second;
755 horizon.push_back(e1);
756 horizon.push_back(e2);
757 horizon.push_back(e3);
763 EdgeIter it = horizon.begin();
764 while (it != horizon.end()) {
765 int occur = count(horizon.begin(), horizon.end(), *it);
767 pair<int, int> removeVal = *it;
768 while (removeVal == *it && !(it == horizon.end()))
769 it = horizon.erase(it);
775 list<myTriangle> newTris;
776 for (EdgeIter it = horizon.begin(); it != horizon.end(); it++) {
777 myTriangle t(it->first, it->second, *farPoint);
778 newTris.push_back(t);
781 vector<myTriangle> trisToProcess;
782 vector<list<int> > newFrontPoints;
783 for (TriIter it = newTris.begin(); it != newTris.end(); it++) {
784 myVec3 t1(xCoords[it->v1], yCoords[it->v1], zCoords[it->v1]);
785 myVec3 t2(xCoords[it->v2], yCoords[it->v2], zCoords[it->v2]);
786 myVec3 t3(xCoords[it->v3], yCoords[it->v3], zCoords[it->v3]);
787 if (isInFront(barycenter, t1, getNorm(t1, t2, t3))) {
796 myVec3 outwardNorm = getNorm(t1, t2, t3);
799 trisToProcess.push_back(tris.back());
801 list<int> newInFront;
802 for (Iter point = pointsInFront.begin(); point != pointsInFront.end();) {
803 myVec3 pointVec(xCoords[*point], yCoords[*point], zCoords[*point]);
804 if (isInFront(pointVec, t1, outwardNorm)) {
805 newInFront.push_back(*point);
806 point = pointsInFront.erase(point);
810 newFrontPoints.push_back(newInFront);
812 vector<myTriangle> allRemoved;
813 for (
int i = 0; i < int(trisToProcess.size()); i++) {
814 if (!newFrontPoints[i].empty()) {
817 if (find(allRemoved.begin(), allRemoved.end(), trisToProcess[i]) == allRemoved.end() && !(trisToProcess[i] == tri)) {
818 vector<myTriangle> removedList = processTriangle(tris, trisToProcess[i], newFrontPoints[i], barycenter, xCoords, yCoords, zCoords);
819 for (
int j = 0; j < int(removedList.size()); j++)
820 allRemoved.push_back(removedList[j]);
827 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
830 "CoarseningVisualization::giftWrap: Gift wrap algorithm input has to have at least 3 points!");
832 #if 1 // TAW's version to determine "minimal" node
834 double min_x = points[0].x;
835 double min_y = points[0].y;
836 for (std::vector<int>::iterator it = nodes.begin(); it != nodes.end(); it++) {
837 int i = it - nodes.begin();
838 if (points[i].x < min_x) min_x = points[i].x;
839 if (points[i].y < min_y) min_y = points[i].y;
844 myVec2 dummy_min(min_x, min_y);
847 std::vector<int> hull;
849 double mindist = distance(min, dummy_min);
850 std::vector<int>::iterator minNode = nodes.begin();
851 for (std::vector<int>::iterator it = nodes.begin(); it != nodes.end(); it++) {
852 int i = it - nodes.begin();
853 if (distance(points[i], dummy_min) < mindist) {
854 mindist = distance(points[i], dummy_min);
859 hull.push_back(*minNode);
860 #else // Brian's code
861 std::vector<int> hull;
862 std::vector<int>::iterator minNode = nodes.begin();
864 for (std::vector<int>::iterator it = nodes.begin(); it != nodes.end(); it++) {
865 int i = it - nodes.begin();
866 if (points[i].x < min.
x || (fabs(points[i].x - min.
x) < 1e-10 && points[i].y < min.
y)) {
871 hull.push_back(*minNode);
874 bool includeMin =
false;
877 std::vector<int>::iterator leftMost = nodes.begin();
878 if (!includeMin && leftMost == minNode) {
881 std::vector<int>::iterator it = leftMost;
883 for (; it != nodes.end(); it++) {
884 if (it == minNode && !includeMin)
886 if (*it == hull.back())
890 myVec2 leftMostVec = points[leftMost - nodes.begin()];
891 myVec2 lastVec(xCoords[hull.back()], yCoords[hull.back()]);
892 myVec2 testNorm = getNorm(vecSubtract(leftMostVec, lastVec));
894 myVec2 itVec(xCoords[*it], yCoords[*it]);
895 double dotProd = dotProduct(testNorm, vecSubtract(itVec, lastVec));
896 if (-1e-8 < dotProd && dotProd < 1e-8) {
899 myVec2 itDist = vecSubtract(itVec, lastVec);
900 myVec2 leftMostDist = vecSubtract(leftMostVec, lastVec);
901 if (fabs(itDist.
x) + fabs(itDist.
y) > fabs(leftMostDist.
x) + fabs(leftMostDist.
y)) {
904 }
else if (dotProd > 0) {
909 if (*leftMost == *minNode)
911 hull.push_back(*leftMost);
919 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
922 vector<int> uniqueNodes = vertices;
923 sort(uniqueNodes.begin(), uniqueNodes.end());
924 vector<int>::iterator newUniqueFineEnd = unique(uniqueNodes.begin(), uniqueNodes.end());
925 uniqueNodes.erase(newUniqueFineEnd, uniqueNodes.end());
928 for (
int i = 0; i < int(vertices.size()); i++) {
930 int hi = uniqueNodes.size() - 1;
932 int search = vertices[i];
934 mid = lo + (hi - lo) / 2;
935 if (uniqueNodes[mid] == search)
937 else if (uniqueNodes[mid] > search)
942 if (uniqueNodes[mid] != search)
943 throw runtime_error(
"Issue in makeUnique_() - a point wasn't found in list.");
949 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
952 const int pos = result.find(replaceWhat);
955 result.replace(pos, replaceWhat.size(), replaceWithWhat);
960 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
962 std::string filenameToWrite = getBaseFileName(numProcs, level, pL);
964 return filenameToWrite;
967 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
969 std::string filenameToWrite = pL.
get<std::string>(
"visualization: output filename");
970 int timeStep = pL.
get<
int>(
"visualization: output file: time step");
971 int iter = pL.
get<
int>(
"visualization: output file: iter");
973 if (filenameToWrite.rfind(
".vtu") == std::string::npos)
974 filenameToWrite.append(
".vtu");
975 if (numProcs > 1 && filenameToWrite.rfind(
"%PROCID") == std::string::npos)
976 filenameToWrite.insert(filenameToWrite.rfind(
".vtu"),
"-proc%PROCID");
981 return filenameToWrite;
984 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
986 std::string filenameToWrite = getBaseFileName(numProcs, level, pL);
987 std::string masterStem = filenameToWrite.substr(0, filenameToWrite.rfind(
".vtu"));
988 masterStem = this->
replaceAll(masterStem,
"%PROCID",
"");
989 std::string pvtuFilename = masterStem +
"-master.pvtu";
993 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
995 std::string styleName =
"PointCloud";
997 std::string indent =
" ";
998 fout <<
"<!--" << styleName <<
" Aggregates Visualization-->" << std::endl;
999 fout <<
"<VTKFile type=\"UnstructuredGrid\" byte_order=\"LittleEndian\">" << std::endl;
1000 fout <<
" <UnstructuredGrid>" << std::endl;
1001 fout <<
" <Piece NumberOfPoints=\"" << uniqueFine.size() <<
"\" NumberOfCells=\"" << geomSizesFine.size() <<
"\">" << std::endl;
1002 fout <<
" <PointData Scalars=\"Node Aggregate Processor\">" << std::endl;
1005 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1007 std::string indent =
" ";
1008 fout <<
" <DataArray type=\"Int32\" Name=\"Node\" format=\"ascii\">" << std::endl;
1011 for (
size_t i = 0; i < uniqueFine.size(); i++) {
1012 if (localIsGlobal) {
1013 fout << uniqueFine[i] <<
" ";
1015 fout << nodeMap->getGlobalElement(uniqueFine[i]) <<
" ";
1021 fout <<
" </DataArray>" << std::endl;
1024 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1026 std::string indent =
" ";
1027 fout <<
" <DataArray type=\"Int32\" Name=\"Aggregate\" format=\"ascii\">" << std::endl;
1029 for (
int i = 0; i < int(uniqueFine.size()); i++) {
1030 fout << myAggOffset + vertex2AggId[uniqueFine[i]] <<
" ";
1036 fout <<
" </DataArray>" << std::endl;
1037 fout <<
" <DataArray type=\"Int32\" Name=\"Processor\" format=\"ascii\">" << std::endl;
1039 for (
int i = 0; i < int(uniqueFine.size()); i++) {
1040 fout << myRank <<
" ";
1046 fout <<
" </DataArray>" << std::endl;
1047 fout <<
" </PointData>" << std::endl;
1050 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1052 std::string indent =
" ";
1053 fout <<
" <Points>" << std::endl;
1054 fout <<
" <DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"ascii\">" << std::endl;
1056 for (
int i = 0; i < int(uniqueFine.size()); i++) {
1057 fout << fx[uniqueFine[i]] <<
" " << fy[uniqueFine[i]] <<
" ";
1061 fout << fz[uniqueFine[i]] <<
" ";
1067 fout <<
" </DataArray>" << std::endl;
1068 fout <<
" </Points>" << std::endl;
1071 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1073 std::string indent =
" ";
1074 fout <<
" <Cells>" << std::endl;
1075 fout <<
" <DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">" << std::endl;
1077 for (
int i = 0; i < int(vertices.size()); i++) {
1078 fout << vertices[i] <<
" ";
1084 fout <<
" </DataArray>" << std::endl;
1085 fout <<
" <DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">" << std::endl;
1088 for (
int i = 0; i < int(geomSize.size()); i++) {
1089 accum += geomSize[i];
1090 fout << accum <<
" ";
1096 fout <<
" </DataArray>" << std::endl;
1097 fout <<
" <DataArray type=\"Int32\" Name=\"types\" format=\"ascii\">" << std::endl;
1099 for (
int i = 0; i < int(geomSize.size()); i++) {
1100 switch (geomSize[i]) {
1118 fout <<
" </DataArray>" << std::endl;
1119 fout <<
" </Cells>" << std::endl;
1122 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1124 fout <<
" </Piece>" << std::endl;
1125 fout <<
" </UnstructuredGrid>" << std::endl;
1126 fout <<
"</VTKFile>" << std::endl;
1129 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1134 pvtu <<
"<VTKFile type=\"PUnstructuredGrid\" byte_order=\"LittleEndian\">" << std::endl;
1135 pvtu <<
" <PUnstructuredGrid GhostLevel=\"0\">" << std::endl;
1136 pvtu <<
" <PPointData Scalars=\"Node Aggregate Processor\">" << std::endl;
1137 pvtu <<
" <PDataArray type=\"Int32\" Name=\"Node\"/>" << std::endl;
1138 pvtu <<
" <PDataArray type=\"Int32\" Name=\"Aggregate\"/>" << std::endl;
1139 pvtu <<
" <PDataArray type=\"Int32\" Name=\"Processor\"/>" << std::endl;
1140 pvtu <<
" </PPointData>" << std::endl;
1141 pvtu <<
" <PPoints>" << std::endl;
1142 pvtu <<
" <PDataArray type=\"Float64\" NumberOfComponents=\"3\"/>" << std::endl;
1143 pvtu <<
" </PPoints>" << std::endl;
1144 for (
int i = 0; i < numProcs; i++) {
1146 pvtu <<
" <Piece Source=\"" <<
replaceAll(baseFname,
"%PROCID",
toString(i)) <<
"\"/>" << std::endl;
1150 for (
int i = 0; i < numProcs; i++) {
1152 pvtu <<
" <Piece Source=\"" << fn.insert(fn.rfind(
".vtu"),
"-finegraph") <<
"\"/>" << std::endl;
1163 pvtu <<
" </PUnstructuredGrid>" << std::endl;
1164 pvtu <<
"</VTKFile>" << std::endl;
1168 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1171 std::ofstream color(
"random-colormap.xml");
1172 color <<
"<ColorMap name=\"MueLu-Random\" space=\"RGB\">" << std::endl;
1175 color <<
" <Point x=\"" << -1 <<
"\" o=\"1\" r=\"1\" g=\"0\" b=\"0\"/>" << std::endl;
1176 color <<
" <Point x=\"" << -2 <<
"\" o=\"1\" r=\"1\" g=\"0.6\" b=\"0\"/>" << std::endl;
1177 color <<
" <Point x=\"" << -3 <<
"\" o=\"1\" r=\"1\" g=\"1\" b=\"0\"/>" << std::endl;
1179 for (
int i = 0; i < 5000; i += 4) {
1180 color <<
" <Point x=\"" << i <<
"\" o=\"1\" r=\"" << (rand() % 50) / 256.0 <<
"\" g=\"" << (rand() % 256) / 256.0 <<
"\" b=\"" << (rand() % 256) / 256.0 <<
"\"/>" << std::endl;
1182 color <<
"</ColorMap>" << std::endl;
1184 }
catch (std::exception& e) {
1186 "VisualizationHelpers::buildColormap: Error while building colormap file: " << e.what());
static void doPointCloud(std::vector< int > &vertices, std::vector< int > &geomSizes, LO numLocalAggs, LO numFineNodes)
MueLu::DefaultLocalOrdinal LocalOrdinal
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
std::string getPVTUFileName(int numProcs, int myRank, int level, const Teuchos::ParameterList &pL) const
void writeFileVTKOpening(std::ofstream &fout, std::vector< int > &uniqueFine, std::vector< int > &geomSizesFine) const
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
std::string replaceAll(std::string result, const std::string &replaceWhat, const std::string &replaceWithWhat) 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)
T & get(const std::string &name, T def_value)
void writeFileVTKData(std::ofstream &fout, std::vector< int > &uniqueFine, LocalOrdinal myAggOffset, ArrayRCP< LocalOrdinal > &vertex2AggId, int myRank) const
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void buildColormap() const
std::vector< int > makeUnique(std::vector< int > &vertices) const
replaces node indices in vertices with compressed unique indices, and returns list of unique points ...
std::string getFileName(int numProcs, int myRank, int level, const Teuchos::ParameterList &pL) const
KOKKOS_INLINE_FUNCTION size_type GetNodeNumVertices() const
Return number of graph vertices.
void writeFileVTKNodes(std::ofstream &fout, std::vector< int > &uniqueFine, Teuchos::RCP< const Map > &nodeMap) const
RCP< ParameterList > GetValidParameterList() 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 sort(View &view, const size_t &size)
static double dotProduct(myVec2 v1, myVec2 v2)
void writePVTU(std::ofstream &pvtu, std::string baseFname, int numProcs, bool bFineEdges=false, bool bCoarseEdges=false) const
std::string getBaseFileName(int numProcs, int level, const Teuchos::ParameterList &pL) const
static void 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)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
MueLu::DefaultGlobalOrdinal GlobalOrdinal
static void doJacks(std::vector< int > &vertices, std::vector< int > &geomSizes, LO numLocalAggs, LO numFineNodes, const std::vector< bool > &isRoot, const ArrayRCP< LO > &vertex2AggId)
static myVec2 getNorm(myVec2 v)
KOKKOS_INLINE_FUNCTION neighbor_vertices_type getNeighborVertices(LO i) const
Return the list of vertices adjacent to the vertex 'v'.
static myVec3 crossProduct(myVec3 v1, myVec3 v2)
static myVec2 vecSubtract(myVec2 v1, myVec2 v2)
static double distance(myVec2 p1, myVec2 p2)
void writeFileVTKCells(std::ofstream &fout, std::vector< int > &uniqueFine, std::vector< LocalOrdinal > &vertices, std::vector< LocalOrdinal > &geomSize) const
static bool isInFront(myVec3 point, myVec3 inPlane, myVec3 n)
Exception throws to report errors in the internal logical of the program.
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 double mymagnitude(myVec2 vec)
static void doGraphEdges(std::vector< int > &vertices, std::vector< int > &geomSizes, Teuchos::RCP< LWGraph > &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)
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)
void writeFileVTKClosing(std::ofstream &fout) const