47 #ifndef MUELU_VISUALIZATIONHELPERS_DEF_HPP_ 48 #define MUELU_VISUALIZATIONHELPERS_DEF_HPP_ 50 #include <Xpetra_Matrix.hpp> 51 #include <Xpetra_CrsMatrixWrap.hpp> 52 #include <Xpetra_ImportFactory.hpp> 53 #include <Xpetra_MultiVectorFactory.hpp> 56 #include "MueLu_Graph.hpp" 65 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
67 RCP<ParameterList> validParamList = rcp(
new ParameterList());
69 validParamList->set< std::string > (
"visualization: output filename",
"viz%LEVELID",
"filename for VTK-style visualization output");
70 validParamList->set<
int > (
"visualization: output file: time step", 0,
"time step variable for output file name");
71 validParamList->set<
int > (
"visualization: output file: iter", 0,
"nonlinear iteration variable for output file name");
72 validParamList->set<std::string> (
"visualization: style",
"Point Cloud",
"style of aggregate visualization for VTK output. Can be 'Point Cloud', 'Jacks', 'Convex Hulls'");
73 validParamList->set<
bool> (
"visualization: build colormap",
false,
"Whether to build a random color map in a separate xml file.");
74 validParamList->set<
bool> (
"visualization: fine graph edges",
false,
"Whether to draw all fine node connections along with the aggregates.");
76 return validParamList;
79 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
81 vertices.reserve(numFineNodes);
82 geomSizes.reserve(numFineNodes);
83 for(LocalOrdinal i = 0; i < numFineNodes; i++)
85 vertices.push_back(i);
86 geomSizes.push_back(1);
90 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
94 vertices.reserve(vertices.size() + 3 * (numFineNodes - numLocalAggs));
95 geomSizes.reserve(vertices.size() + 2 * (numFineNodes - numLocalAggs));
97 for(
int i = 0; i < numLocalAggs; i++)
101 int numInAggFound = 0;
102 for(
int j = 0; j < numFineNodes; j++)
109 if(vertex2AggId[root] == vertex2AggId[j])
111 vertices.push_back(root);
112 vertices.push_back(j);
113 geomSizes.push_back(2);
115 vertices.push_back(j);
116 geomSizes.push_back(1);
126 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
127 void VisualizationHelpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::doConvexHulls2D(std::vector<int>& vertices, std::vector<int>& geomSizes, LO numLocalAggs, LO numFineNodes,
const std::vector<bool>& isRoot,
const ArrayRCP<LO>& vertex2AggId,
const Teuchos::ArrayRCP<const double>& xCoords,
const Teuchos::ArrayRCP<const double>& yCoords,
const Teuchos::ArrayRCP<const double>& zCoords) {
128 for(
int agg = 0; agg < numLocalAggs; agg++) {
129 std::list<int> aggNodes;
130 for(
int i = 0; i < numFineNodes; i++) {
131 if(vertex2AggId[i] == agg)
132 aggNodes.push_back(i);
136 "CoarseningVisualization::doConvexHulls2D: aggregate contains zero nodes!");
137 if(aggNodes.size() == 1) {
138 vertices.push_back(aggNodes.front());
139 geomSizes.push_back(1);
142 if(aggNodes.size() == 2) {
143 vertices.push_back(aggNodes.front());
144 vertices.push_back(aggNodes.back());
145 geomSizes.push_back(2);
149 bool collinear =
true;
151 std::list<int>::iterator it = aggNodes.begin();
152 myVec3 firstPoint(xCoords[*it], yCoords[*it], 0);
154 myVec3 secondPoint(xCoords[*it], yCoords[*it], 0);
156 myVec3 norm1(-(secondPoint.
y - firstPoint.
y), secondPoint.
x - firstPoint.
x, 0);
158 myVec3 thisNorm(yCoords[*it] - firstPoint.
y, firstPoint.
x - xCoords[*it], 0);
160 double temp = thisNorm.
x;
161 thisNorm.
x = thisNorm.
y;
163 double comp = dotProduct(norm1, thisNorm);
164 if(-1e-8 > comp || comp > 1e-8) {
170 while(it != aggNodes.end());
175 std::list<int>::iterator min = aggNodes.begin();
176 std::list<int>::iterator max = aggNodes.begin();
177 for(std::list<int>::iterator it = ++aggNodes.begin(); it != aggNodes.end(); it++) {
178 if(xCoords[*it] < xCoords[*min])
180 else if(xCoords[*it] == xCoords[*min]) {
181 if(yCoords[*it] < yCoords[*min])
184 if(xCoords[*it] > xCoords[*max])
186 else if(xCoords[*it] == xCoords[*max]) {
187 if(yCoords[*it] > yCoords[*max])
192 vertices.push_back(*min);
193 vertices.push_back(*max);
194 geomSizes.push_back(2);
197 std::vector<myVec2> points;
198 std::vector<int> nodes;
199 for(std::list<int>::iterator it = aggNodes.begin(); it != aggNodes.end(); it++) {
200 points.push_back(
myVec2(xCoords[*it], yCoords[*it]));
201 nodes.push_back(*it);
203 std::vector<int> hull = giftWrap(points, nodes, xCoords, yCoords);
204 vertices.reserve(vertices.size() + hull.size());
205 for(
size_t i = 0; i < hull.size(); i++) {
206 vertices.push_back(hull[i]);
208 geomSizes.push_back(hull.size());
212 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
213 void VisualizationHelpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::doConvexHulls3D(std::vector<int>& vertices, std::vector<int>& geomSizes, LO numLocalAggs, LO numFineNodes,
const std::vector<bool>& isRoot,
const ArrayRCP<LO>& vertex2AggId,
const Teuchos::ArrayRCP<const double>& xCoords,
const Teuchos::ArrayRCP<const double>& yCoords,
const Teuchos::ArrayRCP<const double>& zCoords) {
218 typedef std::list<int>::iterator Iter;
219 for(
int agg = 0; agg < numLocalAggs; agg++) {
220 std::list<int> aggNodes;
221 for(
int i = 0; i < numFineNodes; i++) {
222 if(vertex2AggId[i] == agg)
223 aggNodes.push_back(i);
227 "CoarseningVisualization::doConvexHulls3D: aggregate contains zero nodes!");
228 if(aggNodes.size() == 1) {
229 vertices.push_back(aggNodes.front());
230 geomSizes.push_back(1);
232 }
else if(aggNodes.size() == 2) {
233 vertices.push_back(aggNodes.front());
234 vertices.push_back(aggNodes.back());
235 geomSizes.push_back(2);
239 bool areCollinear =
true;
241 Iter it = aggNodes.begin();
242 myVec3 firstVec(xCoords[*it], yCoords[*it], zCoords[*it]);
246 myVec3 secondVec(xCoords[*it], yCoords[*it], zCoords[*it]);
247 comp = vecSubtract(secondVec, firstVec);
250 for(; it != aggNodes.end(); it++) {
251 myVec3 thisVec(xCoords[*it], yCoords[*it], zCoords[*it]);
252 myVec3 cross = crossProduct(vecSubtract(thisVec, firstVec), comp);
253 if(mymagnitude(cross) > 1e-10) {
254 areCollinear =
false;
263 Iter min = aggNodes.begin();
264 Iter max = aggNodes.begin();
265 Iter it = ++aggNodes.begin();
266 for(; it != aggNodes.end(); it++) {
267 if(xCoords[*it] < xCoords[*min]) min = it;
268 else if(xCoords[*it] == xCoords[*min]) {
269 if(yCoords[*it] < yCoords[*min]) min = it;
270 else if(yCoords[*it] == yCoords[*min]) {
271 if(zCoords[*it] < zCoords[*min]) min = it;
274 if(xCoords[*it] > xCoords[*max]) max = it;
275 else if(xCoords[*it] == xCoords[*max]) {
276 if(yCoords[*it] > yCoords[*max]) max = it;
277 else if(yCoords[*it] == yCoords[*max]) {
278 if(zCoords[*it] > zCoords[*max])
283 vertices.push_back(*min);
284 vertices.push_back(*max);
285 geomSizes.push_back(2);
288 bool areCoplanar =
true;
291 Iter vert = aggNodes.begin();
292 myVec3 v1(xCoords[*vert], yCoords[*vert], zCoords[*vert]);
294 myVec3 v2(xCoords[*vert], yCoords[*vert], zCoords[*vert]);
296 myVec3 v3(xCoords[*vert], yCoords[*vert], zCoords[*vert]);
299 while(mymagnitude(crossProduct(vecSubtract(v1, v2), vecSubtract(v1, v3))) < 1e-10) {
301 v3 =
myVec3(xCoords[*vert], yCoords[*vert], zCoords[*vert]);
304 for(; vert != aggNodes.end(); vert++) {
305 myVec3 pt(xCoords[*vert], yCoords[*vert], zCoords[*vert]);
306 if(fabs(pointDistFromTri(pt, v1, v2, v3)) > 1e-12) {
313 myVec3 planeNorm = getNorm(v1, v2, v3);
314 planeNorm.
x = fabs(planeNorm.
x);
315 planeNorm.
y = fabs(planeNorm.
y);
316 planeNorm.
z = fabs(planeNorm.
z);
317 std::vector<myVec2> points;
318 std::vector<int> nodes;
319 if(planeNorm.
x >= planeNorm.
y && planeNorm.
x >= planeNorm.
z) {
321 for(Iter it = aggNodes.begin(); it != aggNodes.end(); it++) {
322 nodes.push_back(*it);
323 points.push_back(
myVec2(yCoords[*it], zCoords[*it]));
326 if(planeNorm.
y >= planeNorm.
x && planeNorm.
y >= planeNorm.
z) {
328 for(Iter it = aggNodes.begin(); it != aggNodes.end(); it++) {
329 nodes.push_back(*it);
330 points.push_back(
myVec2(xCoords[*it], zCoords[*it]));
333 if(planeNorm.
z >= planeNorm.
x && planeNorm.
z >= planeNorm.
y) {
334 for(Iter it = aggNodes.begin(); it != aggNodes.end(); it++) {
335 nodes.push_back(*it);
336 points.push_back(
myVec2(xCoords[*it], yCoords[*it]));
339 std::vector<int> convhull2d = giftWrap(points, nodes, xCoords, yCoords);
340 geomSizes.push_back(convhull2d.size());
341 vertices.reserve(vertices.size() + convhull2d.size());
342 for(
size_t i = 0; i < convhull2d.size(); i++)
343 vertices.push_back(convhull2d[i]);
347 Iter exIt = aggNodes.begin();
348 int extremeSix[] = {*exIt, *exIt, *exIt, *exIt, *exIt, *exIt};
350 for(; exIt != aggNodes.end(); exIt++) {
352 if(xCoords[*it] < xCoords[extremeSix[0]] ||
353 (xCoords[*it] == xCoords[extremeSix[0]] && yCoords[*it] < yCoords[extremeSix[0]]) ||
354 (xCoords[*it] == xCoords[extremeSix[0]] && yCoords[*it] == yCoords[extremeSix[0]] && zCoords[*it] < zCoords[extremeSix[0]]))
356 if(xCoords[*it] > xCoords[extremeSix[1]] ||
357 (xCoords[*it] == xCoords[extremeSix[1]] && yCoords[*it] > yCoords[extremeSix[1]]) ||
358 (xCoords[*it] == xCoords[extremeSix[1]] && yCoords[*it] == yCoords[extremeSix[1]] && zCoords[*it] > zCoords[extremeSix[1]]))
360 if(yCoords[*it] < yCoords[extremeSix[2]] ||
361 (yCoords[*it] == yCoords[extremeSix[2]] && zCoords[*it] < zCoords[extremeSix[2]]) ||
362 (yCoords[*it] == yCoords[extremeSix[2]] && zCoords[*it] == zCoords[extremeSix[2]] && xCoords[*it] < xCoords[extremeSix[2]]))
364 if(yCoords[*it] > yCoords[extremeSix[3]] ||
365 (yCoords[*it] == yCoords[extremeSix[3]] && zCoords[*it] > zCoords[extremeSix[3]]) ||
366 (yCoords[*it] == yCoords[extremeSix[3]] && zCoords[*it] == zCoords[extremeSix[3]] && xCoords[*it] > xCoords[extremeSix[3]]))
368 if(zCoords[*it] < zCoords[extremeSix[4]] ||
369 (zCoords[*it] == zCoords[extremeSix[4]] && xCoords[*it] < xCoords[extremeSix[4]]) ||
370 (zCoords[*it] == zCoords[extremeSix[4]] && xCoords[*it] == xCoords[extremeSix[4]] && yCoords[*it] < yCoords[extremeSix[4]]))
372 if(zCoords[*it] > zCoords[extremeSix[5]] ||
373 (zCoords[*it] == zCoords[extremeSix[5]] && xCoords[*it] > xCoords[extremeSix[5]]) ||
374 (zCoords[*it] == zCoords[extremeSix[5]] && xCoords[*it] == xCoords[extremeSix[5]] && yCoords[*it] > zCoords[extremeSix[5]]))
378 for(
int i = 0; i < 6; i++) {
379 myVec3 thisExtremeVec(xCoords[extremeSix[i]], yCoords[extremeSix[i]], zCoords[extremeSix[i]]);
380 extremeVectors[i] = thisExtremeVec;
385 for(
int i = 0; i < 5; i++) {
386 for(
int j = i + 1; j < 6; j++) {
387 double thisDist = distance(extremeVectors[i], extremeVectors[j]);
388 if(thisDist > maxDist) {
395 std::list<myTriangle> hullBuilding;
397 aggNodes.remove(extremeSix[base1]);
398 aggNodes.remove(extremeSix[base2]);
401 tri.
v1 = extremeSix[base1];
402 tri.
v2 = extremeSix[base2];
406 myVec3 b1 = extremeVectors[base1];
407 myVec3 b2 = extremeVectors[base2];
409 for(Iter node = aggNodes.begin(); node != aggNodes.end(); node++) {
410 myVec3 nodePos(xCoords[*node], yCoords[*node], zCoords[*node]);
411 double dist = mymagnitude(crossProduct(vecSubtract(nodePos, b1), vecSubtract(nodePos, b2))) / mymagnitude(vecSubtract(b2, b1));
419 hullBuilding.push_back(tri);
420 myVec3 b3(xCoords[*thirdNode], yCoords[*thirdNode], zCoords[*thirdNode]);
421 aggNodes.erase(thirdNode);
424 int fourthVertex = -1;
425 for(Iter node = aggNodes.begin(); node != aggNodes.end(); node++) {
426 myVec3 thisNode(xCoords[*node], yCoords[*node], zCoords[*node]);
427 double nodeDist = pointDistFromTri(thisNode, b1, b2, b3);
428 if(nodeDist > maxDist) {
430 fourthVertex = *node;
433 aggNodes.remove(fourthVertex);
434 myVec3 b4(xCoords[fourthVertex], yCoords[fourthVertex], zCoords[fourthVertex]);
437 tri = hullBuilding.front();
438 tri.
v1 = fourthVertex;
439 hullBuilding.push_back(tri);
440 tri = hullBuilding.front();
441 tri.
v2 = fourthVertex;
442 hullBuilding.push_back(tri);
443 tri = hullBuilding.front();
444 tri.
v3 = fourthVertex;
445 hullBuilding.push_back(tri);
447 myVec3 barycenter((b1.
x + b2.
x + b3.
x + b4.
x) / 4.0, (b1.
y + b2.
y + b3.
y + b4.
y) / 4.0, (b1.
z + b2.
z + b3.
z + b4.
z) / 4.0);
448 for(std::list<myTriangle>::iterator tetTri = hullBuilding.begin(); tetTri != hullBuilding.end(); tetTri++) {
449 myVec3 triVert1(xCoords[tetTri->v1], yCoords[tetTri->v1], zCoords[tetTri->v1]);
450 myVec3 triVert2(xCoords[tetTri->v2], yCoords[tetTri->v2], zCoords[tetTri->v2]);
451 myVec3 triVert3(xCoords[tetTri->v3], yCoords[tetTri->v3], zCoords[tetTri->v3]);
452 myVec3 trinorm = getNorm(triVert1, triVert2, triVert3);
453 myVec3 ptInPlane = tetTri == hullBuilding.begin() ? b1 : b4;
454 if(isInFront(barycenter, ptInPlane, trinorm)) {
457 int temp = tetTri->v1;
458 tetTri->v1 = tetTri->v2;
470 for(std::list<myTriangle>::iterator tetTri = hullBuilding.begin(); tetTri != hullBuilding.end(); tetTri++) {
471 myVec3 triVert1(xCoords[tetTri->v1], yCoords[tetTri->v1], zCoords[tetTri->v1]);
472 myVec3 triVert2(xCoords[tetTri->v2], yCoords[tetTri->v2], zCoords[tetTri->v2]);
473 myVec3 triVert3(xCoords[tetTri->v3], yCoords[tetTri->v3], zCoords[tetTri->v3]);
474 trinorms[index] = getNorm(triVert1, triVert2, triVert3);
477 std::list<int> startPoints1;
478 std::list<int> startPoints2;
479 std::list<int> startPoints3;
480 std::list<int> startPoints4;
483 Iter point = aggNodes.begin();
484 while(!aggNodes.empty())
486 myVec3 pointVec(xCoords[*point], yCoords[*point], zCoords[*point]);
489 if(isInFront(pointVec, b1, trinorms[0])) {
490 startPoints1.push_back(*point);
491 point = aggNodes.erase(point);
492 }
else if(isInFront(pointVec, b4, trinorms[1])) {
493 startPoints2.push_back(*point);
494 point = aggNodes.erase(point);
495 }
else if(isInFront(pointVec, b4, trinorms[2])) {
496 startPoints3.push_back(*point);
497 point = aggNodes.erase(point);
498 }
else if(isInFront(pointVec, b4, trinorms[3])) {
499 startPoints4.push_back(*point);
500 point = aggNodes.erase(point);
502 point = aggNodes.erase(point);
507 typedef std::list<myTriangle>::iterator TriIter;
508 TriIter firstTri = hullBuilding.begin();
517 if(!startPoints1.empty())
518 processTriangle(hullBuilding, start1, startPoints1, barycenter, xCoords, yCoords, zCoords);
519 if(!startPoints2.empty())
520 processTriangle(hullBuilding, start2, startPoints2, barycenter, xCoords, yCoords, zCoords);
521 if(!startPoints3.empty())
522 processTriangle(hullBuilding, start3, startPoints3, barycenter, xCoords, yCoords, zCoords);
523 if(!startPoints4.empty())
524 processTriangle(hullBuilding, start4, startPoints4, barycenter, xCoords, yCoords, zCoords);
527 vertices.reserve(vertices.size() + 3 * hullBuilding.size());
528 for(TriIter hullTri = hullBuilding.begin(); hullTri != hullBuilding.end(); hullTri++) {
529 vertices.push_back(hullTri->v1);
530 vertices.push_back(hullTri->v2);
531 vertices.push_back(hullTri->v3);
532 geomSizes.push_back(3);
537 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
539 ArrayView<const Scalar> values;
540 ArrayView<const LocalOrdinal> neighbors;
542 std::vector<std::pair<int, int> > vert1;
544 ArrayView<const LocalOrdinal> indices;
545 for(LocalOrdinal locRow = 0; locRow < LocalOrdinal(G->GetNodeNumVertices()); locRow++) {
546 neighbors = G->getNeighborVertices(locRow);
548 for(
int gEdge = 0; gEdge < int(neighbors.size()); ++gEdge) {
549 vert1.push_back(std::pair<int, int>(locRow, neighbors[gEdge]));
552 for(
size_t i = 0; i < vert1.size(); i ++) {
553 if(vert1[i].first > vert1[i].second) {
554 int temp = vert1[i].first;
555 vert1[i].first = vert1[i].second;
556 vert1[i].second = temp;
559 std::sort(vert1.begin(), vert1.end());
560 std::vector<std::pair<int, int> >::iterator newEnd = unique(vert1.begin(), vert1.end());
561 vert1.erase(newEnd, vert1.end());
563 vertices.reserve(2 * vert1.size());
564 geomSizes.reserve(vert1.size());
565 for(
size_t i = 0; i < vert1.size(); i++) {
566 vertices.push_back(vert1[i].first);
567 vertices.push_back(vert1[i].second);
568 geomSizes.push_back(2);
572 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
575 return myVec3(v1.
y * v2.
z - v1.
z * v2.
y, v1.
z * v2.
x - v1.
x * v2.
z, v1.
x * v2.
y - v1.
y * v2.
x);
578 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
581 return v1.
x * v2.
x + v1.
y * v2.
y;
584 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
587 return v1.
x * v2.
x + v1.
y * v2.
y + v1.
z * v2.
z;
590 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
593 myVec3 rel(point.
x - inPlane.
x, point.
y - inPlane.
y, point.
z - inPlane.
z);
594 return dotProduct(rel, n) > 1e-12 ?
true :
false;
597 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
600 return sqrt(dotProduct(vec, vec));
603 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
606 return sqrt(dotProduct(vec, vec));
609 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
612 return mymagnitude(vecSubtract(p1, p2));
616 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
619 return mymagnitude(vecSubtract(p1, p2));
622 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
628 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
634 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
640 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
643 return crossProduct(vecSubtract(v2, v1), vecSubtract(v3, v1));
647 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
651 myVec3 norm = getNorm(v1, v2, v3);
653 double normScl = mymagnitude(norm);
655 if (normScl > 1e-8) {
659 rv = fabs(dotProduct(norm, vecSubtract(point, v1)));
662 myVec3 test1 = vecSubtract(v3, v1);
663 myVec3 test2 = vecSubtract(v2, v1);
664 bool useTest1 = mymagnitude(test1) > 0.0 ?
true :
false;
665 bool useTest2 = mymagnitude(test2) > 0.0 ?
true :
false;
666 if(useTest1 ==
true) {
667 double normScl1 = mymagnitude(test1);
671 rv = fabs(dotProduct(test1, vecSubtract(point, v1)));
672 }
else if (useTest2 ==
true) {
673 double normScl2 = mymagnitude(test2);
677 rv = fabs(dotProduct(test2, vecSubtract(point, v1)));
680 "VisualizationHelpers::pointDistFromTri: Could not determine the distance of a point to a triangle.");
687 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
693 typedef std::list<int>::iterator Iter;
694 typedef std::list<myTriangle>::iterator TriIter;
695 typedef list<pair<int, int> >::iterator EdgeIter;
698 myVec3 v1(xCoords[tri.
v1], yCoords[tri.
v1], zCoords[tri.
v1]);
699 myVec3 v2(xCoords[tri.
v2], yCoords[tri.
v2], zCoords[tri.
v2]);
700 myVec3 v3(xCoords[tri.
v3], yCoords[tri.
v3], zCoords[tri.
v3]);
702 Iter farPoint = pointsInFront.begin();
703 for(Iter point = pointsInFront.begin(); point != pointsInFront.end(); point++)
705 myVec3 pointVec(xCoords[*point], yCoords[*point], zCoords[*point]);
706 double dist = pointDistFromTri(pointVec, v1, v2, v3);
710 farPointVec = pointVec;
716 vector<myTriangle> visible;
717 for(TriIter it = tris.begin(); it != tris.end();)
719 myVec3 vec1(xCoords[it->v1], yCoords[it->v1], zCoords[it->v1]);
720 myVec3 vec2(xCoords[it->v2], yCoords[it->v2], zCoords[it->v2]);
721 myVec3 vec3(xCoords[it->v3], yCoords[it->v3], zCoords[it->v3]);
722 myVec3 norm = getNorm(vec1, vec2, vec3);
723 if(isInFront(farPointVec, vec1, norm))
725 visible.push_back(*it);
733 list<pair<int, int> > horizon;
736 for(vector<myTriangle>::iterator it = visible.begin(); it != visible.end(); it++)
738 pair<int, int> e1(it->v1, it->v2);
739 pair<int, int> e2(it->v2, it->v3);
740 pair<int, int> e3(it->v1, it->v3);
742 if(e1.first > e1.second)
745 e1.first = e1.second;
748 if(e2.first > e2.second)
751 e2.first = e2.second;
754 if(e3.first > e3.second)
757 e3.first = e3.second;
760 horizon.push_back(e1);
761 horizon.push_back(e2);
762 horizon.push_back(e3);
768 EdgeIter it = horizon.begin();
769 while(it != horizon.end())
771 int occur = count(horizon.begin(), horizon.end(), *it);
774 pair<int, int> removeVal = *it;
775 while(removeVal == *it && !(it == horizon.end()))
776 it = horizon.erase(it);
783 list<myTriangle> newTris;
784 for(EdgeIter it = horizon.begin(); it != horizon.end(); it++)
786 myTriangle t(it->first, it->second, *farPoint);
787 newTris.push_back(t);
790 vector<myTriangle> trisToProcess;
791 vector<list<int> > newFrontPoints;
792 for(TriIter it = newTris.begin(); it != newTris.end(); it++)
794 myVec3 t1(xCoords[it->v1], yCoords[it->v1], zCoords[it->v1]);
795 myVec3 t2(xCoords[it->v2], yCoords[it->v2], zCoords[it->v2]);
796 myVec3 t3(xCoords[it->v3], yCoords[it->v3], zCoords[it->v3]);
797 if(isInFront(barycenter, t1, getNorm(t1, t2, t3)))
807 myVec3 outwardNorm = getNorm(t1, t2, t3);
810 trisToProcess.push_back(tris.back());
812 list<int> newInFront;
813 for(Iter point = pointsInFront.begin(); point != pointsInFront.end();)
815 myVec3 pointVec(xCoords[*point], yCoords[*point], zCoords[*point]);
816 if(isInFront(pointVec, t1, outwardNorm))
818 newInFront.push_back(*point);
819 point = pointsInFront.erase(point);
824 newFrontPoints.push_back(newInFront);
826 vector<myTriangle> allRemoved;
827 for(
int i = 0; i < int(trisToProcess.size()); i++)
829 if(!newFrontPoints[i].empty())
833 if(find(allRemoved.begin(), allRemoved.end(), trisToProcess[i]) == allRemoved.end() && !(trisToProcess[i] == tri))
835 vector<myTriangle> removedList = processTriangle(tris, trisToProcess[i], newFrontPoints[i], barycenter, xCoords, yCoords, zCoords);
836 for(
int j = 0; j < int(removedList.size()); j++)
837 allRemoved.push_back(removedList[j]);
844 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
847 "CoarseningVisualization::giftWrap: Gift wrap algorithm input has to have at least 3 points!");
849 #if 1 // TAW's version to determine "minimal" node 851 double min_x =points[0].x;
852 double min_y =points[0].y;
853 for(std::vector<int>::iterator it = nodes.begin(); it != nodes.end(); it++) {
854 int i = it - nodes.begin();
855 if(points[i].x < min_x) min_x = points[i].x;
856 if(points[i].y < min_y) min_y = points[i].y;
861 myVec2 dummy_min(min_x, min_y);
864 std::vector<int> hull;
866 double mindist = distance(min,dummy_min);
867 std::vector<int>::iterator minNode = nodes.begin();
868 for(std::vector<int>::iterator it = nodes.begin(); it != nodes.end(); it++) {
869 int i = it - nodes.begin();
870 if(distance(points[i],dummy_min) < mindist) {
871 mindist = distance(points[i],dummy_min);
876 hull.push_back(*minNode);
877 #else // Brian's code 878 std::vector<int> hull;
879 std::vector<int>::iterator minNode = nodes.begin();
881 for(std::vector<int>::iterator it = nodes.begin(); it != nodes.end(); it++)
883 int i = it - nodes.begin();
884 if(points[i].x < min.
x || (fabs(points[i].x - min.
x) < 1e-10 && points[i].y < min.
y))
890 hull.push_back(*minNode);
893 bool includeMin =
false;
897 std::vector<int>::iterator leftMost = nodes.begin();
898 if(!includeMin && leftMost == minNode)
902 std::vector<int>::iterator it = leftMost;
904 for(; it != nodes.end(); it++)
906 if(it == minNode && !includeMin)
908 if(*it == hull.back())
912 myVec2 leftMostVec = points[leftMost - nodes.begin()];
913 myVec2 lastVec(xCoords[hull.back()], yCoords[hull.back()]);
914 myVec2 testNorm = getNorm(vecSubtract(leftMostVec, lastVec));
916 myVec2 itVec(xCoords[*it], yCoords[*it]);
917 double dotProd = dotProduct(testNorm, vecSubtract(itVec, lastVec));
918 if(-1e-8 < dotProd && dotProd < 1e-8)
922 myVec2 itDist = vecSubtract(itVec, lastVec);
923 myVec2 leftMostDist = vecSubtract(leftMostVec, lastVec);
924 if(fabs(itDist.
x) + fabs(itDist.
y) > fabs(leftMostDist.
x) + fabs(leftMostDist.
y)) {
928 else if(dotProd > 0) {
934 if(*leftMost == *minNode)
936 hull.push_back(*leftMost);
944 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
948 vector<int> uniqueNodes = vertices;
949 sort(uniqueNodes.begin(), uniqueNodes.end());
950 vector<int>::iterator newUniqueFineEnd = unique(uniqueNodes.begin(), uniqueNodes.end());
951 uniqueNodes.erase(newUniqueFineEnd, uniqueNodes.end());
954 for(
int i = 0; i < int(vertices.size()); i++)
957 int hi = uniqueNodes.size() - 1;
959 int search = vertices[i];
962 mid = lo + (hi - lo) / 2;
963 if(uniqueNodes[mid] == search)
965 else if(uniqueNodes[mid] > search)
970 if(uniqueNodes[mid] != search)
971 throw runtime_error(
"Issue in makeUnique_() - a point wasn't found in list.");
977 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
980 const int pos = result.find(replaceWhat);
983 result.replace(pos, replaceWhat.size(), replaceWithWhat);
988 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
990 std::string filenameToWrite = getBaseFileName(numProcs, level, pL);
992 return filenameToWrite;
995 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
997 std::string filenameToWrite = pL.get<std::string>(
"visualization: output filename");
998 int timeStep = pL.get<
int>(
"visualization: output file: time step");
999 int iter = pL.get<
int>(
"visualization: output file: iter");
1001 if(filenameToWrite.rfind(
".vtu") == std::string::npos)
1002 filenameToWrite.append(
".vtu");
1003 if(numProcs > 1 && filenameToWrite.rfind(
"%PROCID") == std::string::npos)
1004 filenameToWrite.insert(filenameToWrite.rfind(
".vtu"),
"-proc%PROCID");
1007 filenameToWrite = this->
replaceAll(filenameToWrite,
"%TIMESTEP",
toString(timeStep));
1009 return filenameToWrite;
1012 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1014 std::string filenameToWrite = getBaseFileName(numProcs, level, pL);
1015 std::string masterStem = filenameToWrite.substr(0, filenameToWrite.rfind(
".vtu"));
1016 masterStem = this->
replaceAll(masterStem,
"%PROCID",
"");
1017 std::string pvtuFilename = masterStem +
"-master.pvtu";
1018 return pvtuFilename;
1021 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1023 std::string styleName =
"PointCloud";
1025 std::string indent =
" ";
1026 fout <<
"<!--" << styleName <<
" Aggregates Visualization-->" << std::endl;
1027 fout <<
"<VTKFile type=\"UnstructuredGrid\" byte_order=\"LittleEndian\">" << std::endl;
1028 fout <<
" <UnstructuredGrid>" << std::endl;
1029 fout <<
" <Piece NumberOfPoints=\"" << uniqueFine.size() <<
"\" NumberOfCells=\"" << geomSizesFine.size() <<
"\">" << std::endl;
1030 fout <<
" <PointData Scalars=\"Node Aggregate Processor\">" << std::endl;
1033 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1035 std::string indent =
" ";
1036 fout <<
" <DataArray type=\"Int32\" Name=\"Node\" format=\"ascii\">" << std::endl;
1038 bool localIsGlobal = GlobalOrdinal(nodeMap->getGlobalNumElements()) == GlobalOrdinal(nodeMap->getNodeNumElements());
1039 for(
size_t i = 0; i < uniqueFine.size(); i++)
1043 fout << uniqueFine[i] <<
" ";
1046 fout << nodeMap->getGlobalElement(uniqueFine[i]) <<
" ";
1048 fout << std::endl << indent;
1051 fout <<
" </DataArray>" << std::endl;
1054 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1056 std::string indent =
" ";
1057 fout <<
" <DataArray type=\"Int32\" Name=\"Aggregate\" format=\"ascii\">" << std::endl;
1059 for(
int i = 0; i < int(uniqueFine.size()); i++)
1061 fout << myAggOffset + vertex2AggId[uniqueFine[i]] <<
" ";
1063 fout << std::endl << indent;
1066 fout <<
" </DataArray>" << std::endl;
1067 fout <<
" <DataArray type=\"Int32\" Name=\"Processor\" format=\"ascii\">" << std::endl;
1069 for(
int i = 0; i < int(uniqueFine.size()); i++)
1071 fout << myRank <<
" ";
1073 fout << std::endl << indent;
1076 fout <<
" </DataArray>" << std::endl;
1077 fout <<
" </PointData>" << std::endl;
1080 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1082 std::string indent =
" ";
1083 fout <<
" <Points>" << std::endl;
1084 fout <<
" <DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"ascii\">" << std::endl;
1086 for(
int i = 0; i < int(uniqueFine.size()); i++)
1088 fout << fx[uniqueFine[i]] <<
" " << fy[uniqueFine[i]] <<
" ";
1092 fout << fz[uniqueFine[i]] <<
" ";
1094 fout << std::endl << indent;
1097 fout <<
" </DataArray>" << std::endl;
1098 fout <<
" </Points>" << std::endl;
1101 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1103 std::string indent =
" ";
1104 fout <<
" <Cells>" << std::endl;
1105 fout <<
" <DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">" << std::endl;
1107 for(
int i = 0; i < int(vertices.size()); i++)
1109 fout << vertices[i] <<
" ";
1111 fout << std::endl << indent;
1114 fout <<
" </DataArray>" << std::endl;
1115 fout <<
" <DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">" << std::endl;
1118 for(
int i = 0; i < int(geomSize.size()); i++)
1120 accum += geomSize[i];
1121 fout << accum <<
" ";
1123 fout << std::endl << indent;
1126 fout <<
" </DataArray>" << std::endl;
1127 fout <<
" <DataArray type=\"Int32\" Name=\"types\" format=\"ascii\">" << std::endl;
1129 for(
int i = 0; i < int(geomSize.size()); i++)
1146 fout << std::endl << indent;
1149 fout <<
" </DataArray>" << std::endl;
1150 fout <<
" </Cells>" << std::endl;
1153 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1155 fout <<
" </Piece>" << std::endl;
1156 fout <<
" </UnstructuredGrid>" << std::endl;
1157 fout <<
"</VTKFile>" << std::endl;
1161 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1166 pvtu <<
"<VTKFile type=\"PUnstructuredGrid\" byte_order=\"LittleEndian\">" << std::endl;
1167 pvtu <<
" <PUnstructuredGrid GhostLevel=\"0\">" << std::endl;
1168 pvtu <<
" <PPointData Scalars=\"Node Aggregate Processor\">" << std::endl;
1169 pvtu <<
" <PDataArray type=\"Int32\" Name=\"Node\"/>" << std::endl;
1170 pvtu <<
" <PDataArray type=\"Int32\" Name=\"Aggregate\"/>" << std::endl;
1171 pvtu <<
" <PDataArray type=\"Int32\" Name=\"Processor\"/>" << std::endl;
1172 pvtu <<
" </PPointData>" << std::endl;
1173 pvtu <<
" <PPoints>" << std::endl;
1174 pvtu <<
" <PDataArray type=\"Float64\" NumberOfComponents=\"3\"/>" << std::endl;
1175 pvtu <<
" </PPoints>" << std::endl;
1176 for(
int i = 0; i < numProcs; i++) {
1178 pvtu <<
" <Piece Source=\"" <<
replaceAll(baseFname,
"%PROCID",
toString(i)) <<
"\"/>" << std::endl;
1183 for(
int i = 0; i < numProcs; i++)
1186 pvtu <<
" <Piece Source=\"" << fn.insert(fn.rfind(
".vtu"),
"-finegraph") <<
"\"/>" << std::endl;
1197 pvtu <<
" </PUnstructuredGrid>" << std::endl;
1198 pvtu <<
"</VTKFile>" << std::endl;
1202 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1205 std::ofstream color(
"random-colormap.xml");
1206 color <<
"<ColorMap name=\"MueLu-Random\" space=\"RGB\">" << std::endl;
1209 color <<
" <Point x=\"" << -1 <<
"\" o=\"1\" r=\"1\" g=\"0\" b=\"0\"/>" << std::endl;
1210 color <<
" <Point x=\"" << -2 <<
"\" o=\"1\" r=\"1\" g=\"0.6\" b=\"0\"/>" << std::endl;
1211 color <<
" <Point x=\"" << -3 <<
"\" o=\"1\" r=\"1\" g=\"1\" b=\"0\"/>" << std::endl;
1213 for(
int i = 0; i < 5000; i += 4) {
1214 color <<
" <Point x=\"" << i <<
"\" o=\"1\" r=\"" << (rand() % 50) / 256.0 <<
"\" g=\"" << (rand() % 256) / 256.0 <<
"\" b=\"" << (rand() % 256) / 256.0 <<
"\"/>" << std::endl;
1216 color <<
"</ColorMap>" << std::endl;
1219 catch(std::exception& e) {
1221 "VisualizationHelpers::buildColormap: Error while building colormap file: " << e.what());
static void doPointCloud(std::vector< int > &vertices, std::vector< int > &geomSizes, LO numLocalAggs, LO numFineNodes)
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
std::string getPVTUFileName(int numProcs, int myRank, int level, const Teuchos::ParameterList &pL) const
void writeFileVTKOpening(std::ofstream &fout, std::vector< int > &uniqueFine, std::vector< int > &geomSizesFine) const
std::string replaceAll(std::string result, const std::string &replaceWhat, const std::string &replaceWithWhat) const
void writeFileVTKData(std::ofstream &fout, std::vector< int > &uniqueFine, LocalOrdinal myAggOffset, ArrayRCP< LocalOrdinal > &vertex2AggId, int myRank) const
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
Namespace for MueLu classes and methods.
void writeFileVTKNodes(std::ofstream &fout, std::vector< int > &uniqueFine, Teuchos::RCP< const Map > &nodeMap) const
RCP< ParameterList > GetValidParameterList() const
static double pointDistFromTri(myVec3 point, myVec3 v1, myVec3 v2, myVec3 v3)
static double dotProduct(myVec2 v1, myVec2 v2)
void writePVTU(std::ofstream &pvtu, std::string baseFname, int numProcs, bool bFineEdges=false, bool bCoarseEdges=false) const
std::string getBaseFileName(int numProcs, int level, const Teuchos::ParameterList &pL) const
static void doGraphEdges(std::vector< int > &vertices, std::vector< int > &geomSizes, Teuchos::RCP< GraphBase > &G, Teuchos::ArrayRCP< const double > &fx, Teuchos::ArrayRCP< const double > &fy, Teuchos::ArrayRCP< const double > &fz)
static std::vector< int > giftWrap(std::vector< myVec2 > &points, std::vector< int > &nodes, const Teuchos::ArrayRCP< const double > &xCoords, const Teuchos::ArrayRCP< const double > &yCoords)
static void doJacks(std::vector< int > &vertices, std::vector< int > &geomSizes, LO numLocalAggs, LO numFineNodes, const std::vector< bool > &isRoot, const ArrayRCP< LO > &vertex2AggId)
static myVec2 getNorm(myVec2 v)
static void doConvexHulls2D(std::vector< int > &vertices, std::vector< int > &geomSizes, LO numLocalAggs, LO numFineNodes, const std::vector< bool > &isRoot, const ArrayRCP< LO > &vertex2AggId, const Teuchos::ArrayRCP< const double > &xCoords, const Teuchos::ArrayRCP< const double > &yCoords, const Teuchos::ArrayRCP< const double > &zCoords)
static myVec3 crossProduct(myVec3 v1, myVec3 v2)
void writeFileVTKCoordinates(std::ofstream &fout, std::vector< int > &uniqueFine, Teuchos::ArrayRCP< const double > &fx, Teuchos::ArrayRCP< const double > &fy, Teuchos::ArrayRCP< const double > &fz, int dim) const
static myVec2 vecSubtract(myVec2 v1, myVec2 v2)
static std::vector< myTriangle > processTriangle(std::list< myTriangle > &tris, myTriangle tri, std::list< int > &pointsInFront, myVec3 &barycenter, const Teuchos::ArrayRCP< const double > &xCoords, const Teuchos::ArrayRCP< const double > &yCoords, const Teuchos::ArrayRCP< const double > &zCoords)
static void doConvexHulls3D(std::vector< int > &vertices, std::vector< int > &geomSizes, LO numLocalAggs, LO numFineNodes, const std::vector< bool > &isRoot, const ArrayRCP< LO > &vertex2AggId, const Teuchos::ArrayRCP< const double > &xCoords, const Teuchos::ArrayRCP< const double > &yCoords, const Teuchos::ArrayRCP< const double > &zCoords)
static double distance(myVec2 p1, myVec2 p2)
void writeFileVTKCells(std::ofstream &fout, std::vector< int > &uniqueFine, std::vector< LocalOrdinal > &vertices, std::vector< LocalOrdinal > &geomSize) const
static bool isInFront(myVec3 point, myVec3 inPlane, myVec3 n)
Exception throws to report errors in the internal logical of the program.
static double mymagnitude(myVec2 vec)
void replaceAll(std::string &str, const std::string &from, const std::string &to)
void writeFileVTKClosing(std::ofstream &fout) const