47 #ifndef MUELU_VISUALIZATIONHELPERS_DEF_HPP_
48 #define MUELU_VISUALIZATIONHELPERS_DEF_HPP_
56 #include "MueLu_Graph.hpp"
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>
73 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
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");
79 validParamList->
set<
int > (
"visualization: output file: iter", 0,
"nonlinear iteration variable for output file name");
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.");
84 return validParamList;
87 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
89 vertices.reserve(numFineNodes);
90 geomSizes.reserve(numFineNodes);
91 for(LocalOrdinal i = 0; i < numFineNodes; i++)
93 vertices.push_back(i);
94 geomSizes.push_back(1);
98 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
102 vertices.reserve(vertices.size() + 3 * (numFineNodes - numLocalAggs));
103 geomSizes.reserve(vertices.size() + 2 * (numFineNodes - numLocalAggs));
105 for(
int i = 0; i < numLocalAggs; i++)
109 int numInAggFound = 0;
110 for(
int j = 0; j < numFineNodes; j++)
117 if(vertex2AggId[root] == vertex2AggId[j])
119 vertices.push_back(root);
120 vertices.push_back(j);
121 geomSizes.push_back(2);
123 vertices.push_back(j);
124 geomSizes.push_back(1);
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>& ,
const ArrayRCP<LO>& vertex2AggId,
const Teuchos::ArrayRCP<
const typename Teuchos::ScalarTraits<Scalar>::coordinateType>& xCoords,
const Teuchos::ArrayRCP<
const typename Teuchos::ScalarTraits<Scalar>::coordinateType>& yCoords) {
136 for(
int agg = 0; agg < numLocalAggs; agg++) {
137 std::list<int> aggNodes;
138 for(
int i = 0; i < numFineNodes; i++) {
139 if(vertex2AggId[i] == agg)
140 aggNodes.push_back(i);
144 "CoarseningVisualization::doConvexHulls2D: aggregate contains zero nodes!");
145 if(aggNodes.size() == 1) {
146 vertices.push_back(aggNodes.front());
147 geomSizes.push_back(1);
150 if(aggNodes.size() == 2) {
151 vertices.push_back(aggNodes.front());
152 vertices.push_back(aggNodes.back());
153 geomSizes.push_back(2);
157 bool collinear =
true;
159 std::list<int>::iterator it = aggNodes.begin();
160 myVec3 firstPoint(xCoords[*it], yCoords[*it], 0);
162 myVec3 secondPoint(xCoords[*it], yCoords[*it], 0);
164 myVec3 norm1(-(secondPoint.
y - firstPoint.
y), secondPoint.
x - firstPoint.
x, 0);
166 myVec3 thisNorm(yCoords[*it] - firstPoint.
y, firstPoint.
x - xCoords[*it], 0);
168 double temp = thisNorm.
x;
169 thisNorm.
x = thisNorm.
y;
171 double comp = dotProduct(norm1, thisNorm);
172 if(-1e-8 > comp || comp > 1e-8) {
178 while(it != aggNodes.end());
183 std::list<int>::iterator min = aggNodes.begin();
184 std::list<int>::iterator max = aggNodes.begin();
185 for(std::list<int>::iterator it = ++aggNodes.begin(); it != aggNodes.end(); it++) {
186 if(xCoords[*it] < xCoords[*min])
188 else if(xCoords[*it] == xCoords[*min]) {
189 if(yCoords[*it] < yCoords[*min])
192 if(xCoords[*it] > xCoords[*max])
194 else if(xCoords[*it] == xCoords[*max]) {
195 if(yCoords[*it] > yCoords[*max])
200 vertices.push_back(*min);
201 vertices.push_back(*max);
202 geomSizes.push_back(2);
205 std::vector<myVec2> points;
206 std::vector<int> nodes;
207 for(std::list<int>::iterator it = aggNodes.begin(); it != aggNodes.end(); it++) {
208 points.push_back(
myVec2(xCoords[*it], yCoords[*it]));
209 nodes.push_back(*it);
211 std::vector<int> hull = giftWrap(points, nodes, xCoords, yCoords);
212 vertices.reserve(vertices.size() + hull.size());
213 for(
size_t i = 0; i < hull.size(); i++) {
214 vertices.push_back(hull[i]);
216 geomSizes.push_back(hull.size());
220 #ifdef HAVE_MUELU_CGAL
221 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
222 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) {
224 typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
225 typedef K::Point_2 Point_2;
227 for(
int agg = 0; agg < numLocalAggs; agg++) {
228 std::vector<int> aggNodes;
229 std::vector<Point_2> aggPoints;
230 for(
int i = 0; i < numFineNodes; i++) {
231 if(vertex2AggId[i] == agg) {
232 Point_2 p(xCoords[i], yCoords[i]);
233 aggPoints.push_back(p);
234 aggNodes.push_back(i);
239 "CoarseningVisualization::doCGALConvexHulls2D: aggregate contains zero nodes!");
240 if(aggNodes.size() == 1) {
241 vertices.push_back(aggNodes.front());
242 geomSizes.push_back(1);
245 if(aggNodes.size() == 2) {
246 vertices.push_back(aggNodes.front());
247 vertices.push_back(aggNodes.back());
248 geomSizes.push_back(2);
252 bool collinear =
true;
254 std::vector<int>::iterator it = aggNodes.begin();
255 myVec3 firstPoint(xCoords[*it], yCoords[*it], 0);
257 myVec3 secondPoint(xCoords[*it], yCoords[*it], 0);
259 myVec3 norm1(-(secondPoint.
y - firstPoint.
y), secondPoint.
x - firstPoint.
x, 0);
261 myVec3 thisNorm(yCoords[*it] - firstPoint.
y, firstPoint.
x - xCoords[*it], 0);
263 double temp = thisNorm.
x;
264 thisNorm.
x = thisNorm.
y;
266 double comp = dotProduct(norm1, thisNorm);
267 if(-1e-8 > comp || comp > 1e-8) {
273 while(it != aggNodes.end());
278 std::vector<int>::iterator min = aggNodes.begin();
279 std::vector<int>::iterator max = aggNodes.begin();
280 for(std::vector<int>::iterator it = ++aggNodes.begin(); it != aggNodes.end(); it++) {
281 if(xCoords[*it] < xCoords[*min])
283 else if(xCoords[*it] == xCoords[*min]) {
284 if(yCoords[*it] < yCoords[*min])
287 if(xCoords[*it] > xCoords[*max])
289 else if(xCoords[*it] == xCoords[*max]) {
290 if(yCoords[*it] > yCoords[*max])
295 vertices.push_back(*min);
296 vertices.push_back(*max);
297 geomSizes.push_back(2);
302 std::vector<Point_2> result;
303 CGAL::convex_hull_2( aggPoints.begin(), aggPoints.end(), std::back_inserter(result) );
305 for (
size_t r = 0; r < result.size(); r++) {
307 for(
size_t l = 0; l < aggNodes.size(); l++)
309 if(fabs(result[r].x() - xCoords[aggNodes[l]]) < 1e-12 &&
310 fabs(result[r].y() - yCoords[aggNodes[l]]) < 1e-12)
312 vertices.push_back(aggNodes[l]);
318 geomSizes.push_back(result.size());
325 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
326 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) {
331 typedef std::list<int>::iterator Iter;
332 for(
int agg = 0; agg < numLocalAggs; agg++) {
333 std::list<int> aggNodes;
334 for(
int i = 0; i < numFineNodes; i++) {
335 if(vertex2AggId[i] == agg)
336 aggNodes.push_back(i);
340 "CoarseningVisualization::doConvexHulls3D: aggregate contains zero nodes!");
341 if(aggNodes.size() == 1) {
342 vertices.push_back(aggNodes.front());
343 geomSizes.push_back(1);
345 }
else if(aggNodes.size() == 2) {
346 vertices.push_back(aggNodes.front());
347 vertices.push_back(aggNodes.back());
348 geomSizes.push_back(2);
352 bool areCollinear =
true;
354 Iter it = aggNodes.begin();
355 myVec3 firstVec(xCoords[*it], yCoords[*it], zCoords[*it]);
359 myVec3 secondVec(xCoords[*it], yCoords[*it], zCoords[*it]);
360 comp = vecSubtract(secondVec, firstVec);
363 for(; it != aggNodes.end(); it++) {
364 myVec3 thisVec(xCoords[*it], yCoords[*it], zCoords[*it]);
365 myVec3 cross = crossProduct(vecSubtract(thisVec, firstVec), comp);
366 if(mymagnitude(cross) > 1e-10) {
367 areCollinear =
false;
376 Iter min = aggNodes.begin();
377 Iter max = aggNodes.begin();
378 Iter it = ++aggNodes.begin();
379 for(; it != aggNodes.end(); it++) {
380 if(xCoords[*it] < xCoords[*min]) min = it;
381 else if(xCoords[*it] == xCoords[*min]) {
382 if(yCoords[*it] < yCoords[*min]) min = it;
383 else if(yCoords[*it] == yCoords[*min]) {
384 if(zCoords[*it] < zCoords[*min]) min = it;
387 if(xCoords[*it] > xCoords[*max]) max = it;
388 else if(xCoords[*it] == xCoords[*max]) {
389 if(yCoords[*it] > yCoords[*max]) max = it;
390 else if(yCoords[*it] == yCoords[*max]) {
391 if(zCoords[*it] > zCoords[*max])
396 vertices.push_back(*min);
397 vertices.push_back(*max);
398 geomSizes.push_back(2);
401 bool areCoplanar =
true;
404 Iter vert = aggNodes.begin();
405 myVec3 v1(xCoords[*vert], yCoords[*vert], zCoords[*vert]);
407 myVec3 v2(xCoords[*vert], yCoords[*vert], zCoords[*vert]);
409 myVec3 v3(xCoords[*vert], yCoords[*vert], zCoords[*vert]);
412 while(mymagnitude(crossProduct(vecSubtract(v1, v2), vecSubtract(v1, v3))) < 1e-10) {
414 v3 =
myVec3(xCoords[*vert], yCoords[*vert], zCoords[*vert]);
417 for(; vert != aggNodes.end(); vert++) {
418 myVec3 pt(xCoords[*vert], yCoords[*vert], zCoords[*vert]);
419 if(fabs(pointDistFromTri(pt, v1, v2, v3)) > 1e-12) {
426 myVec3 planeNorm = getNorm(v1, v2, v3);
427 planeNorm.
x = fabs(planeNorm.
x);
428 planeNorm.
y = fabs(planeNorm.
y);
429 planeNorm.
z = fabs(planeNorm.
z);
430 std::vector<myVec2> points;
431 std::vector<int> nodes;
432 if(planeNorm.
x >= planeNorm.
y && planeNorm.
x >= planeNorm.
z) {
434 for(Iter it = aggNodes.begin(); it != aggNodes.end(); it++) {
435 nodes.push_back(*it);
436 points.push_back(
myVec2(yCoords[*it], zCoords[*it]));
439 if(planeNorm.
y >= planeNorm.
x && planeNorm.
y >= planeNorm.
z) {
441 for(Iter it = aggNodes.begin(); it != aggNodes.end(); it++) {
442 nodes.push_back(*it);
443 points.push_back(
myVec2(xCoords[*it], zCoords[*it]));
446 if(planeNorm.
z >= planeNorm.
x && planeNorm.
z >= planeNorm.
y) {
447 for(Iter it = aggNodes.begin(); it != aggNodes.end(); it++) {
448 nodes.push_back(*it);
449 points.push_back(
myVec2(xCoords[*it], yCoords[*it]));
452 std::vector<int> convhull2d = giftWrap(points, nodes, xCoords, yCoords);
453 geomSizes.push_back(convhull2d.size());
454 vertices.reserve(vertices.size() + convhull2d.size());
455 for(
size_t i = 0; i < convhull2d.size(); i++)
456 vertices.push_back(convhull2d[i]);
460 Iter exIt = aggNodes.begin();
461 int extremeSix[] = {*exIt, *exIt, *exIt, *exIt, *exIt, *exIt};
463 for(; exIt != aggNodes.end(); exIt++) {
465 if(xCoords[*it] < xCoords[extremeSix[0]] ||
466 (xCoords[*it] == xCoords[extremeSix[0]] && yCoords[*it] < yCoords[extremeSix[0]]) ||
467 (xCoords[*it] == xCoords[extremeSix[0]] && yCoords[*it] == yCoords[extremeSix[0]] && zCoords[*it] < zCoords[extremeSix[0]]))
469 if(xCoords[*it] > xCoords[extremeSix[1]] ||
470 (xCoords[*it] == xCoords[extremeSix[1]] && yCoords[*it] > yCoords[extremeSix[1]]) ||
471 (xCoords[*it] == xCoords[extremeSix[1]] && yCoords[*it] == yCoords[extremeSix[1]] && zCoords[*it] > zCoords[extremeSix[1]]))
473 if(yCoords[*it] < yCoords[extremeSix[2]] ||
474 (yCoords[*it] == yCoords[extremeSix[2]] && zCoords[*it] < zCoords[extremeSix[2]]) ||
475 (yCoords[*it] == yCoords[extremeSix[2]] && zCoords[*it] == zCoords[extremeSix[2]] && xCoords[*it] < xCoords[extremeSix[2]]))
477 if(yCoords[*it] > yCoords[extremeSix[3]] ||
478 (yCoords[*it] == yCoords[extremeSix[3]] && zCoords[*it] > zCoords[extremeSix[3]]) ||
479 (yCoords[*it] == yCoords[extremeSix[3]] && zCoords[*it] == zCoords[extremeSix[3]] && xCoords[*it] > xCoords[extremeSix[3]]))
481 if(zCoords[*it] < zCoords[extremeSix[4]] ||
482 (zCoords[*it] == zCoords[extremeSix[4]] && xCoords[*it] < xCoords[extremeSix[4]]) ||
483 (zCoords[*it] == zCoords[extremeSix[4]] && xCoords[*it] == xCoords[extremeSix[4]] && yCoords[*it] < yCoords[extremeSix[4]]))
485 if(zCoords[*it] > zCoords[extremeSix[5]] ||
486 (zCoords[*it] == zCoords[extremeSix[5]] && xCoords[*it] > xCoords[extremeSix[5]]) ||
487 (zCoords[*it] == zCoords[extremeSix[5]] && xCoords[*it] == xCoords[extremeSix[5]] && yCoords[*it] > zCoords[extremeSix[5]]))
491 for(
int i = 0; i < 6; i++) {
492 myVec3 thisExtremeVec(xCoords[extremeSix[i]], yCoords[extremeSix[i]], zCoords[extremeSix[i]]);
493 extremeVectors[i] = thisExtremeVec;
498 for(
int i = 0; i < 5; i++) {
499 for(
int j = i + 1; j < 6; j++) {
500 double thisDist = distance(extremeVectors[i], extremeVectors[j]);
501 if(thisDist > maxDist) {
508 std::list<myTriangle> hullBuilding;
510 aggNodes.remove(extremeSix[base1]);
511 aggNodes.remove(extremeSix[base2]);
514 tri.
v1 = extremeSix[base1];
515 tri.
v2 = extremeSix[base2];
519 myVec3 b1 = extremeVectors[base1];
520 myVec3 b2 = extremeVectors[base2];
522 for(Iter node = aggNodes.begin(); node != aggNodes.end(); node++) {
523 myVec3 nodePos(xCoords[*node], yCoords[*node], zCoords[*node]);
524 double dist = mymagnitude(crossProduct(vecSubtract(nodePos, b1), vecSubtract(nodePos, b2))) / mymagnitude(vecSubtract(b2, b1));
532 hullBuilding.push_back(tri);
533 myVec3 b3(xCoords[*thirdNode], yCoords[*thirdNode], zCoords[*thirdNode]);
534 aggNodes.erase(thirdNode);
537 int fourthVertex = -1;
538 for(Iter node = aggNodes.begin(); node != aggNodes.end(); node++) {
539 myVec3 thisNode(xCoords[*node], yCoords[*node], zCoords[*node]);
540 double nodeDist = pointDistFromTri(thisNode, b1, b2, b3);
541 if(nodeDist > maxDist) {
543 fourthVertex = *node;
546 aggNodes.remove(fourthVertex);
547 myVec3 b4(xCoords[fourthVertex], yCoords[fourthVertex], zCoords[fourthVertex]);
550 tri = hullBuilding.front();
551 tri.
v1 = fourthVertex;
552 hullBuilding.push_back(tri);
553 tri = hullBuilding.front();
554 tri.
v2 = fourthVertex;
555 hullBuilding.push_back(tri);
556 tri = hullBuilding.front();
557 tri.
v3 = fourthVertex;
558 hullBuilding.push_back(tri);
560 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);
561 for(std::list<myTriangle>::iterator tetTri = hullBuilding.begin(); tetTri != hullBuilding.end(); tetTri++) {
562 myVec3 triVert1(xCoords[tetTri->v1], yCoords[tetTri->v1], zCoords[tetTri->v1]);
563 myVec3 triVert2(xCoords[tetTri->v2], yCoords[tetTri->v2], zCoords[tetTri->v2]);
564 myVec3 triVert3(xCoords[tetTri->v3], yCoords[tetTri->v3], zCoords[tetTri->v3]);
565 myVec3 trinorm = getNorm(triVert1, triVert2, triVert3);
566 myVec3 ptInPlane = tetTri == hullBuilding.begin() ? b1 : b4;
567 if(isInFront(barycenter, ptInPlane, trinorm)) {
570 int temp = tetTri->v1;
571 tetTri->v1 = tetTri->v2;
583 for(std::list<myTriangle>::iterator tetTri = hullBuilding.begin(); tetTri != hullBuilding.end(); tetTri++) {
584 myVec3 triVert1(xCoords[tetTri->v1], yCoords[tetTri->v1], zCoords[tetTri->v1]);
585 myVec3 triVert2(xCoords[tetTri->v2], yCoords[tetTri->v2], zCoords[tetTri->v2]);
586 myVec3 triVert3(xCoords[tetTri->v3], yCoords[tetTri->v3], zCoords[tetTri->v3]);
587 trinorms[index] = getNorm(triVert1, triVert2, triVert3);
590 std::list<int> startPoints1;
591 std::list<int> startPoints2;
592 std::list<int> startPoints3;
593 std::list<int> startPoints4;
596 Iter point = aggNodes.begin();
597 while(!aggNodes.empty())
599 myVec3 pointVec(xCoords[*point], yCoords[*point], zCoords[*point]);
602 if(isInFront(pointVec, b1, trinorms[0])) {
603 startPoints1.push_back(*point);
604 point = aggNodes.erase(point);
605 }
else if(isInFront(pointVec, b4, trinorms[1])) {
606 startPoints2.push_back(*point);
607 point = aggNodes.erase(point);
608 }
else if(isInFront(pointVec, b4, trinorms[2])) {
609 startPoints3.push_back(*point);
610 point = aggNodes.erase(point);
611 }
else if(isInFront(pointVec, b4, trinorms[3])) {
612 startPoints4.push_back(*point);
613 point = aggNodes.erase(point);
615 point = aggNodes.erase(point);
620 typedef std::list<myTriangle>::iterator TriIter;
621 TriIter firstTri = hullBuilding.begin();
630 if(!startPoints1.empty())
631 processTriangle(hullBuilding, start1, startPoints1, barycenter, xCoords, yCoords, zCoords);
632 if(!startPoints2.empty())
633 processTriangle(hullBuilding, start2, startPoints2, barycenter, xCoords, yCoords, zCoords);
634 if(!startPoints3.empty())
635 processTriangle(hullBuilding, start3, startPoints3, barycenter, xCoords, yCoords, zCoords);
636 if(!startPoints4.empty())
637 processTriangle(hullBuilding, start4, startPoints4, barycenter, xCoords, yCoords, zCoords);
640 vertices.reserve(vertices.size() + 3 * hullBuilding.size());
641 for(TriIter hullTri = hullBuilding.begin(); hullTri != hullBuilding.end(); hullTri++) {
642 vertices.push_back(hullTri->v1);
643 vertices.push_back(hullTri->v2);
644 vertices.push_back(hullTri->v3);
645 geomSizes.push_back(3);
650 #ifdef HAVE_MUELU_CGAL
651 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
652 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) {
654 typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
655 typedef K::Point_3 Point_3;
656 typedef CGAL::Polyhedron_3<K> Polyhedron_3;
657 typedef std::vector<int>::iterator Iter;
658 for(
int agg = 0; agg < numLocalAggs; agg++) {
659 std::vector<int> aggNodes;
660 std::vector<Point_3> aggPoints;
661 for(
int i = 0; i < numFineNodes; i++) {
662 if(vertex2AggId[i] == agg) {
663 Point_3 p(xCoords[i], yCoords[i], zCoords[i]);
664 aggPoints.push_back(p);
665 aggNodes.push_back(i);
670 "CoarseningVisualization::doCGALConvexHulls3D: aggregate contains zero nodes!");
671 if(aggNodes.size() == 1) {
672 vertices.push_back(aggNodes.front());
673 geomSizes.push_back(1);
675 }
else if(aggNodes.size() == 2) {
676 vertices.push_back(aggNodes.front());
677 vertices.push_back(aggNodes.back());
678 geomSizes.push_back(2);
682 bool areCollinear =
true;
684 Iter it = aggNodes.begin();
685 myVec3 firstVec(xCoords[*it], yCoords[*it], zCoords[*it]);
689 myVec3 secondVec(xCoords[*it], yCoords[*it], zCoords[*it]);
690 comp = vecSubtract(secondVec, firstVec);
693 for(; it != aggNodes.end(); it++) {
694 myVec3 thisVec(xCoords[*it], yCoords[*it], zCoords[*it]);
695 myVec3 cross = crossProduct(vecSubtract(thisVec, firstVec), comp);
696 if(mymagnitude(cross) > 1e-8) {
697 areCollinear =
false;
706 Iter min = aggNodes.begin();
707 Iter max = aggNodes.begin();
708 Iter it = ++aggNodes.begin();
709 for(; it != aggNodes.end(); it++) {
710 if(xCoords[*it] < xCoords[*min]) min = it;
711 else if(xCoords[*it] == xCoords[*min]) {
712 if(yCoords[*it] < yCoords[*min]) min = it;
713 else if(yCoords[*it] == yCoords[*min]) {
714 if(zCoords[*it] < zCoords[*min]) min = it;
717 if(xCoords[*it] > xCoords[*max]) max = it;
718 else if(xCoords[*it] == xCoords[*max]) {
719 if(yCoords[*it] > yCoords[*max]) max = it;
720 else if(yCoords[*it] == yCoords[*max]) {
721 if(zCoords[*it] > zCoords[*max])
726 vertices.push_back(*min);
727 vertices.push_back(*max);
728 geomSizes.push_back(2);
734 CGAL::convex_hull_3( aggPoints.begin(), aggPoints.end(), result);
737 Polyhedron_3::Facet_iterator fi;
738 for (fi = result.facets_begin(); fi != result.facets_end(); fi++) {
739 int cntVertInAgg = 0;
740 Polyhedron_3::Halfedge_around_facet_const_circulator hit = fi->facet_begin();
742 const Point_3 & pp = hit->vertex()->point();
744 for(
size_t l = 0; l < aggNodes.size(); l++)
746 if(fabs(pp.x() - xCoords[aggNodes[l]]) < 1e-12 &&
747 fabs(pp.y() - yCoords[aggNodes[l]]) < 1e-12 &&
748 fabs(pp.z() - zCoords[aggNodes[l]]) < 1e-12)
750 vertices.push_back(aggNodes[l]);
755 }
while (++hit != fi->facet_begin());
756 geomSizes.push_back(cntVertInAgg);
764 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
769 std::vector<std::pair<int, int> > vert1;
772 for(LocalOrdinal locRow = 0; locRow < LocalOrdinal(G->
GetNodeNumVertices()); locRow++) {
775 for(
int gEdge = 0; gEdge < int(neighbors.
size()); ++gEdge) {
776 vert1.push_back(std::pair<int, int>(locRow, neighbors[gEdge]));
779 for(
size_t i = 0; i < vert1.size(); i ++) {
780 if(vert1[i].first > vert1[i].second) {
781 int temp = vert1[i].first;
782 vert1[i].first = vert1[i].second;
783 vert1[i].second = temp;
786 std::sort(vert1.begin(), vert1.end());
787 std::vector<std::pair<int, int> >::iterator newEnd = unique(vert1.begin(), vert1.end());
788 vert1.erase(newEnd, vert1.end());
790 vertices.reserve(2 * vert1.size());
791 geomSizes.reserve(vert1.size());
792 for(
size_t i = 0; i < vert1.size(); i++) {
793 vertices.push_back(vert1[i].first);
794 vertices.push_back(vert1[i].second);
795 geomSizes.push_back(2);
799 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
802 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);
805 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
808 return v1.
x * v2.
x + v1.
y * v2.
y;
811 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
814 return v1.
x * v2.
x + v1.
y * v2.
y + v1.
z * v2.
z;
817 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
820 myVec3 rel(point.
x - inPlane.
x, point.
y - inPlane.
y, point.
z - inPlane.
z);
821 return dotProduct(rel, n) > 1e-12 ?
true :
false;
824 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
827 return sqrt(dotProduct(vec, vec));
830 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
833 return sqrt(dotProduct(vec, vec));
836 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
839 return mymagnitude(vecSubtract(p1, p2));
843 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
846 return mymagnitude(vecSubtract(p1, p2));
849 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
855 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
861 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
867 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
870 return crossProduct(vecSubtract(v2, v1), vecSubtract(v3, v1));
874 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
878 myVec3 norm = getNorm(v1, v2, v3);
880 double normScl = mymagnitude(norm);
882 if (normScl > 1e-8) {
886 rv = fabs(dotProduct(norm, vecSubtract(point, v1)));
889 myVec3 test1 = vecSubtract(v3, v1);
890 myVec3 test2 = vecSubtract(v2, v1);
891 bool useTest1 = mymagnitude(test1) > 0.0 ?
true :
false;
892 bool useTest2 = mymagnitude(test2) > 0.0 ?
true :
false;
893 if(useTest1 ==
true) {
894 double normScl1 = mymagnitude(test1);
898 rv = fabs(dotProduct(test1, vecSubtract(point, v1)));
899 }
else if (useTest2 ==
true) {
900 double normScl2 = mymagnitude(test2);
904 rv = fabs(dotProduct(test2, vecSubtract(point, v1)));
907 "VisualizationHelpers::pointDistFromTri: Could not determine the distance of a point to a triangle.");
914 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
920 typedef std::list<int>::iterator Iter;
921 typedef std::list<myTriangle>::iterator TriIter;
922 typedef list<pair<int, int> >::iterator EdgeIter;
925 myVec3 v1(xCoords[tri.
v1], yCoords[tri.
v1], zCoords[tri.
v1]);
926 myVec3 v2(xCoords[tri.
v2], yCoords[tri.
v2], zCoords[tri.
v2]);
927 myVec3 v3(xCoords[tri.
v3], yCoords[tri.
v3], zCoords[tri.
v3]);
929 Iter farPoint = pointsInFront.begin();
930 for(Iter point = pointsInFront.begin(); point != pointsInFront.end(); point++)
932 myVec3 pointVec(xCoords[*point], yCoords[*point], zCoords[*point]);
933 double dist = pointDistFromTri(pointVec, v1, v2, v3);
937 farPointVec = pointVec;
943 vector<myTriangle> visible;
944 for(TriIter it = tris.begin(); it != tris.end();)
946 myVec3 vec1(xCoords[it->v1], yCoords[it->v1], zCoords[it->v1]);
947 myVec3 vec2(xCoords[it->v2], yCoords[it->v2], zCoords[it->v2]);
948 myVec3 vec3(xCoords[it->v3], yCoords[it->v3], zCoords[it->v3]);
949 myVec3 norm = getNorm(vec1, vec2, vec3);
950 if(isInFront(farPointVec, vec1, norm))
952 visible.push_back(*it);
960 list<pair<int, int> > horizon;
963 for(vector<myTriangle>::iterator it = visible.begin(); it != visible.end(); it++)
965 pair<int, int> e1(it->v1, it->v2);
966 pair<int, int> e2(it->v2, it->v3);
967 pair<int, int> e3(it->v1, it->v3);
969 if(e1.first > e1.second)
972 e1.first = e1.second;
975 if(e2.first > e2.second)
978 e2.first = e2.second;
981 if(e3.first > e3.second)
984 e3.first = e3.second;
987 horizon.push_back(e1);
988 horizon.push_back(e2);
989 horizon.push_back(e3);
995 EdgeIter it = horizon.begin();
996 while(it != horizon.end())
998 int occur = count(horizon.begin(), horizon.end(), *it);
1001 pair<int, int> removeVal = *it;
1002 while(removeVal == *it && !(it == horizon.end()))
1003 it = horizon.erase(it);
1010 list<myTriangle> newTris;
1011 for(EdgeIter it = horizon.begin(); it != horizon.end(); it++)
1013 myTriangle t(it->first, it->second, *farPoint);
1014 newTris.push_back(t);
1017 vector<myTriangle> trisToProcess;
1018 vector<list<int> > newFrontPoints;
1019 for(TriIter it = newTris.begin(); it != newTris.end(); it++)
1021 myVec3 t1(xCoords[it->v1], yCoords[it->v1], zCoords[it->v1]);
1022 myVec3 t2(xCoords[it->v2], yCoords[it->v2], zCoords[it->v2]);
1023 myVec3 t3(xCoords[it->v3], yCoords[it->v3], zCoords[it->v3]);
1024 if(isInFront(barycenter, t1, getNorm(t1, t2, t3)))
1034 myVec3 outwardNorm = getNorm(t1, t2, t3);
1036 tris.push_back(*it);
1037 trisToProcess.push_back(tris.back());
1039 list<int> newInFront;
1040 for(Iter point = pointsInFront.begin(); point != pointsInFront.end();)
1042 myVec3 pointVec(xCoords[*point], yCoords[*point], zCoords[*point]);
1043 if(isInFront(pointVec, t1, outwardNorm))
1045 newInFront.push_back(*point);
1046 point = pointsInFront.erase(point);
1051 newFrontPoints.push_back(newInFront);
1053 vector<myTriangle> allRemoved;
1054 for(
int i = 0; i < int(trisToProcess.size()); i++)
1056 if(!newFrontPoints[i].empty())
1060 if(find(allRemoved.begin(), allRemoved.end(), trisToProcess[i]) == allRemoved.end() && !(trisToProcess[i] == tri))
1062 vector<myTriangle> removedList = processTriangle(tris, trisToProcess[i], newFrontPoints[i], barycenter, xCoords, yCoords, zCoords);
1063 for(
int j = 0; j < int(removedList.size()); j++)
1064 allRemoved.push_back(removedList[j]);
1071 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1074 "CoarseningVisualization::giftWrap: Gift wrap algorithm input has to have at least 3 points!");
1076 #if 1 // TAW's version to determine "minimal" node
1078 double min_x =points[0].x;
1079 double min_y =points[0].y;
1080 for(std::vector<int>::iterator it = nodes.begin(); it != nodes.end(); it++) {
1081 int i = it - nodes.begin();
1082 if(points[i].x < min_x) min_x = points[i].x;
1083 if(points[i].y < min_y) min_y = points[i].y;
1088 myVec2 dummy_min(min_x, min_y);
1091 std::vector<int> hull;
1093 double mindist = distance(min,dummy_min);
1094 std::vector<int>::iterator minNode = nodes.begin();
1095 for(std::vector<int>::iterator it = nodes.begin(); it != nodes.end(); it++) {
1096 int i = it - nodes.begin();
1097 if(distance(points[i],dummy_min) < mindist) {
1098 mindist = distance(points[i],dummy_min);
1103 hull.push_back(*minNode);
1104 #else // Brian's code
1105 std::vector<int> hull;
1106 std::vector<int>::iterator minNode = nodes.begin();
1108 for(std::vector<int>::iterator it = nodes.begin(); it != nodes.end(); it++)
1110 int i = it - nodes.begin();
1111 if(points[i].x < min.
x || (fabs(points[i].x - min.
x) < 1e-10 && points[i].y < min.
y))
1117 hull.push_back(*minNode);
1120 bool includeMin =
false;
1124 std::vector<int>::iterator leftMost = nodes.begin();
1125 if(!includeMin && leftMost == minNode)
1129 std::vector<int>::iterator it = leftMost;
1131 for(; it != nodes.end(); it++)
1133 if(it == minNode && !includeMin)
1135 if(*it == hull.back())
1139 myVec2 leftMostVec = points[leftMost - nodes.begin()];
1140 myVec2 lastVec(xCoords[hull.back()], yCoords[hull.back()]);
1141 myVec2 testNorm = getNorm(vecSubtract(leftMostVec, lastVec));
1143 myVec2 itVec(xCoords[*it], yCoords[*it]);
1144 double dotProd = dotProduct(testNorm, vecSubtract(itVec, lastVec));
1145 if(-1e-8 < dotProd && dotProd < 1e-8)
1149 myVec2 itDist = vecSubtract(itVec, lastVec);
1150 myVec2 leftMostDist = vecSubtract(leftMostVec, lastVec);
1151 if(fabs(itDist.
x) + fabs(itDist.
y) > fabs(leftMostDist.
x) + fabs(leftMostDist.
y)) {
1155 else if(dotProd > 0) {
1161 if(*leftMost == *minNode)
1163 hull.push_back(*leftMost);
1171 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1174 using namespace std;
1175 vector<int> uniqueNodes = vertices;
1176 sort(uniqueNodes.begin(), uniqueNodes.end());
1177 vector<int>::iterator newUniqueFineEnd = unique(uniqueNodes.begin(), uniqueNodes.end());
1178 uniqueNodes.erase(newUniqueFineEnd, uniqueNodes.end());
1181 for(
int i = 0; i < int(vertices.size()); i++)
1184 int hi = uniqueNodes.size() - 1;
1186 int search = vertices[i];
1189 mid = lo + (hi - lo) / 2;
1190 if(uniqueNodes[mid] == search)
1192 else if(uniqueNodes[mid] > search)
1197 if(uniqueNodes[mid] != search)
1198 throw runtime_error(
"Issue in makeUnique_() - a point wasn't found in list.");
1204 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1207 const int pos = result.find(replaceWhat);
1210 result.replace(pos, replaceWhat.size(), replaceWithWhat);
1215 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1217 std::string filenameToWrite = getBaseFileName(numProcs, level, pL);
1219 return filenameToWrite;
1222 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1224 std::string filenameToWrite = pL.
get<std::string>(
"visualization: output filename");
1225 int timeStep = pL.
get<
int>(
"visualization: output file: time step");
1226 int iter = pL.
get<
int>(
"visualization: output file: iter");
1228 if(filenameToWrite.rfind(
".vtu") == std::string::npos)
1229 filenameToWrite.append(
".vtu");
1230 if(numProcs > 1 && filenameToWrite.rfind(
"%PROCID") == std::string::npos)
1231 filenameToWrite.insert(filenameToWrite.rfind(
".vtu"),
"-proc%PROCID");
1234 filenameToWrite = this->
replaceAll(filenameToWrite,
"%TIMESTEP",
toString(timeStep));
1236 return filenameToWrite;
1239 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1241 std::string filenameToWrite = getBaseFileName(numProcs, level, pL);
1242 std::string masterStem = filenameToWrite.substr(0, filenameToWrite.rfind(
".vtu"));
1243 masterStem = this->
replaceAll(masterStem,
"%PROCID",
"");
1244 std::string pvtuFilename = masterStem +
"-master.pvtu";
1245 return pvtuFilename;
1248 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1250 std::string styleName =
"PointCloud";
1252 std::string indent =
" ";
1253 fout <<
"<!--" << styleName <<
" Aggregates Visualization-->" << std::endl;
1254 fout <<
"<VTKFile type=\"UnstructuredGrid\" byte_order=\"LittleEndian\">" << std::endl;
1255 fout <<
" <UnstructuredGrid>" << std::endl;
1256 fout <<
" <Piece NumberOfPoints=\"" << uniqueFine.size() <<
"\" NumberOfCells=\"" << geomSizesFine.size() <<
"\">" << std::endl;
1257 fout <<
" <PointData Scalars=\"Node Aggregate Processor\">" << std::endl;
1260 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1262 std::string indent =
" ";
1263 fout <<
" <DataArray type=\"Int32\" Name=\"Node\" format=\"ascii\">" << std::endl;
1265 bool localIsGlobal = GlobalOrdinal(nodeMap->getGlobalNumElements()) == GlobalOrdinal(nodeMap->getNodeNumElements());
1266 for(
size_t i = 0; i < uniqueFine.size(); i++)
1270 fout << uniqueFine[i] <<
" ";
1273 fout << nodeMap->getGlobalElement(uniqueFine[i]) <<
" ";
1275 fout << std::endl << indent;
1278 fout <<
" </DataArray>" << std::endl;
1281 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1283 std::string indent =
" ";
1284 fout <<
" <DataArray type=\"Int32\" Name=\"Aggregate\" format=\"ascii\">" << std::endl;
1286 for(
int i = 0; i < int(uniqueFine.size()); i++)
1288 fout << myAggOffset + vertex2AggId[uniqueFine[i]] <<
" ";
1290 fout << std::endl << indent;
1293 fout <<
" </DataArray>" << std::endl;
1294 fout <<
" <DataArray type=\"Int32\" Name=\"Processor\" format=\"ascii\">" << std::endl;
1296 for(
int i = 0; i < int(uniqueFine.size()); i++)
1298 fout << myRank <<
" ";
1300 fout << std::endl << indent;
1303 fout <<
" </DataArray>" << std::endl;
1304 fout <<
" </PointData>" << std::endl;
1307 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1309 std::string indent =
" ";
1310 fout <<
" <Points>" << std::endl;
1311 fout <<
" <DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"ascii\">" << std::endl;
1313 for(
int i = 0; i < int(uniqueFine.size()); i++)
1315 fout << fx[uniqueFine[i]] <<
" " << fy[uniqueFine[i]] <<
" ";
1319 fout << fz[uniqueFine[i]] <<
" ";
1321 fout << std::endl << indent;
1324 fout <<
" </DataArray>" << std::endl;
1325 fout <<
" </Points>" << std::endl;
1328 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1330 std::string indent =
" ";
1331 fout <<
" <Cells>" << std::endl;
1332 fout <<
" <DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">" << std::endl;
1334 for(
int i = 0; i < int(vertices.size()); i++)
1336 fout << vertices[i] <<
" ";
1338 fout << std::endl << indent;
1341 fout <<
" </DataArray>" << std::endl;
1342 fout <<
" <DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">" << std::endl;
1345 for(
int i = 0; i < int(geomSize.size()); i++)
1347 accum += geomSize[i];
1348 fout << accum <<
" ";
1350 fout << std::endl << indent;
1353 fout <<
" </DataArray>" << std::endl;
1354 fout <<
" <DataArray type=\"Int32\" Name=\"types\" format=\"ascii\">" << std::endl;
1356 for(
int i = 0; i < int(geomSize.size()); i++)
1373 fout << std::endl << indent;
1376 fout <<
" </DataArray>" << std::endl;
1377 fout <<
" </Cells>" << std::endl;
1380 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1382 fout <<
" </Piece>" << std::endl;
1383 fout <<
" </UnstructuredGrid>" << std::endl;
1384 fout <<
"</VTKFile>" << std::endl;
1388 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1393 pvtu <<
"<VTKFile type=\"PUnstructuredGrid\" byte_order=\"LittleEndian\">" << std::endl;
1394 pvtu <<
" <PUnstructuredGrid GhostLevel=\"0\">" << std::endl;
1395 pvtu <<
" <PPointData Scalars=\"Node Aggregate Processor\">" << std::endl;
1396 pvtu <<
" <PDataArray type=\"Int32\" Name=\"Node\"/>" << std::endl;
1397 pvtu <<
" <PDataArray type=\"Int32\" Name=\"Aggregate\"/>" << std::endl;
1398 pvtu <<
" <PDataArray type=\"Int32\" Name=\"Processor\"/>" << std::endl;
1399 pvtu <<
" </PPointData>" << std::endl;
1400 pvtu <<
" <PPoints>" << std::endl;
1401 pvtu <<
" <PDataArray type=\"Float64\" NumberOfComponents=\"3\"/>" << std::endl;
1402 pvtu <<
" </PPoints>" << std::endl;
1403 for(
int i = 0; i < numProcs; i++) {
1405 pvtu <<
" <Piece Source=\"" <<
replaceAll(baseFname,
"%PROCID",
toString(i)) <<
"\"/>" << std::endl;
1410 for(
int i = 0; i < numProcs; i++)
1413 pvtu <<
" <Piece Source=\"" << fn.insert(fn.rfind(
".vtu"),
"-finegraph") <<
"\"/>" << std::endl;
1424 pvtu <<
" </PUnstructuredGrid>" << std::endl;
1425 pvtu <<
"</VTKFile>" << std::endl;
1429 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1432 std::ofstream color(
"random-colormap.xml");
1433 color <<
"<ColorMap name=\"MueLu-Random\" space=\"RGB\">" << std::endl;
1436 color <<
" <Point x=\"" << -1 <<
"\" o=\"1\" r=\"1\" g=\"0\" b=\"0\"/>" << std::endl;
1437 color <<
" <Point x=\"" << -2 <<
"\" o=\"1\" r=\"1\" g=\"0.6\" b=\"0\"/>" << std::endl;
1438 color <<
" <Point x=\"" << -3 <<
"\" o=\"1\" r=\"1\" g=\"1\" b=\"0\"/>" << std::endl;
1440 for(
int i = 0; i < 5000; i += 4) {
1441 color <<
" <Point x=\"" << i <<
"\" o=\"1\" r=\"" << (rand() % 50) / 256.0 <<
"\" g=\"" << (rand() % 256) / 256.0 <<
"\" b=\"" << (rand() % 256) / 256.0 <<
"\"/>" << std::endl;
1443 color <<
"</ColorMap>" << std::endl;
1446 catch(std::exception& e) {
1448 "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
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)
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)
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
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)
virtual size_t GetNodeNumVertices() const =0
Return number of vertices owned by the calling node.
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
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 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)
static void doGraphEdges(std::vector< int > &vertices, std::vector< int > &geomSizes, Teuchos::RCP< GraphBase > &G, Teuchos::ArrayRCP< const typename Teuchos::ScalarTraits< Scalar >::coordinateType > &fx, Teuchos::ArrayRCP< const typename Teuchos::ScalarTraits< Scalar >::coordinateType > &fy, Teuchos::ArrayRCP< const typename Teuchos::ScalarTraits< Scalar >::coordinateType > &fz)
static void doJacks(std::vector< int > &vertices, std::vector< int > &geomSizes, LO numLocalAggs, LO numFineNodes, const std::vector< bool > &isRoot, const ArrayRCP< LO > &vertex2AggId)
static myVec2 getNorm(myVec2 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)
void replaceAll(std::string &str, const std::string &from, const std::string &to)
virtual Teuchos::ArrayView< const LocalOrdinal > getNeighborVertices(LocalOrdinal v) const =0
Return the list of vertices adjacent to the vertex 'v'.
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