Sundance 2.1

00001 /*
00002  * SundanceSurf3DCalc.cpp
00003  *
00004  *  Created on: Oct 25, 2011
00005  *      Author: benk
00006  */
00007 
00008 #include "SundanceSurf3DCalc.hpp"
00009 
00010 using namespace Sundance;
00011 
00012 const int SundanceSurf3DCalc::edegIndex[12][2] = { {0,1} , {0,2} , {0,4} , {1,3} , {1,5} , {2,3} , {2,6} , {3,7} ,
00013                                {4,5} , {4,6} , {5,7} , {6,7} };
00014                                               //       0           1           2           3            4             5
00015 const int SundanceSurf3DCalc::faceEdges[6][4]  = { {0,1,3,5} , {0,2,4,8} , {1,2,6,9} , {3,4,7,10} , {5,6,7,11} , {8,9,10,11} };
00016 
00017                                               //     0       1       2       3       4       5       6       7      8       9       10      11
00018 const int SundanceSurf3DCalc::edgeFaces[12][2] = { {0,1} , {0,2} , {1,2} , {0,3} , {1,3} , {0,4} , {2,4} , {3,4} , {1,5} , {2,5} , {3,5} , {4,5} };
00019 
00020 //                                                            0   1   2   3   4   5   6   7   8   9   10  11
00021 const int SundanceSurf3DCalc::edgeNeighboredges[12][12] = { { 1 , 1 , 1 , 1 , 1 , 1 , 0 , 0 , 1 , 0 , 0 , 0} ,   /* 0 */
00022                                                         { 1 , 1 , 1 , 1 , 0 , 1 , 1 , 0 , 0 , 1 , 0 , 0} ,   /* 1 */
00023                                                         { 1 , 1 , 1 , 0 , 1 , 0 , 1 , 0 , 1 , 1 , 0 , 0} ,   /* 2 */
00024                                                         { 1 , 1 , 0 , 1 , 1 , 1 , 0 , 1 , 0 , 0 , 1 , 0} ,   /* 3 */
00025                                                         { 1 , 0 , 1 , 1 , 1 , 0 , 0 , 1 , 1 , 0 , 1 , 0} ,   /* 4 */
00026                                                         { 1 , 1 , 0 , 1 , 0 , 1 , 1 , 1 , 0 , 0 , 0 , 1} ,   /* 5 */
00027                                                         { 0 , 1 , 1 , 0 , 0 , 1 , 1 , 1 , 0 , 1 , 0 , 1} ,   /* 6 */
00028                                                         { 0 , 0 , 0 , 1 , 1 , 1 , 1 , 1 , 0 , 0 , 1 , 1} ,   /* 7 */
00029                                                         { 1 , 0 , 1 , 0 , 1 , 0 , 0 , 0 , 1 , 1 , 1 , 1} ,   /* 8 */
00030                                                         { 0 , 1 , 1 , 0 , 0 , 0 , 1 , 0 , 1 , 1 , 1 , 1} ,   /* 9 */
00031                                                         { 0 , 0 , 0 , 1 , 1 , 0 , 0 , 1 , 1 , 1 , 1 , 1} ,   /* 10 */
00032                                                         { 0 , 0 , 0 , 0 , 0 , 1 , 1 , 1 , 1 , 1 , 1 , 1} , };/* 11 */
00033 
00034 void SundanceSurf3DCalc::getCurveQuadPoints(CellType  maxCellType ,
00035                                int maxCellLID ,
00036                                const Mesh& mesh ,
00037                    const ParametrizedCurve& paramCurve,
00038                    const Array<Point>& brickPoints,
00039                    Array<Point>& intersectPoints ,
00040                    Array<int>& edgeIndex,
00041                    Array<int>& triangleIndex){
00042   TEUCHOS_TEST_FOR_EXCEPTION( maxCellType != BrickCell ,
00043           std::runtime_error , " maxCellType must be a BrickCell: " << BrickCell << " but it is:" << maxCellType);
00044 
00045   int nrPoints , edge , t_inters_points = 0;
00046   int verb = 2;
00047   int tmp_i;
00048   int faceIntPoint[6] = { 0 , 0 , 0 , 0 , 0 , 0 };
00049   int edgePointI[12] = { -1 , -1 , -1 , -1 , -1 , -1 , -1 , -1 , -1 , -1 , -1 , -1 };
00050 
00051 // -------------------- get the intersection point according to the policy -------------------
00052 
00053   // get the intersection points between the edges of the brick and the surface
00054   for (edge = 0 ; edge < 12 ; edge++){
00055       Array<Point> result(0);
00056     paramCurve.returnIntersectPoints(brickPoints[edegIndex[edge][0]], brickPoints[edegIndex[edge][1]], nrPoints , result);
00057     // test if we have intersection point, and take this only if there is one
00058       if (nrPoints == 1){
00059         // we test if the impacted faces have already two intersection points, if yes then, ignore the intersection point
00060         // ignore this intersection point
00061         if ( (faceIntPoint[ edgeFaces[edge][0]] < 2) &&  (faceIntPoint[ edgeFaces[edge][1] ] < 2 )) {
00062            faceIntPoint[ edgeFaces[edge][0] ]++; faceIntPoint[ edgeFaces[edge][1] ]++;
00063            intersectPoints.resize( t_inters_points + 1 );
00064            edgeIndex.resize( t_inters_points + 1 );
00065            intersectPoints[t_inters_points] = result[0];
00066            edgeIndex[t_inters_points] = edge;
00067            edgePointI[edge] = t_inters_points;
00068            //SUNDANCE_MSG3(verb, "found Int. point [" << t_inters_points <<"]=" << result[0]);
00069            t_inters_points++;
00070         } else {
00071           SUNDANCE_MSG1(verb, "WARNING One face had already two intersection points :"
00072               << faceIntPoint[ edgeFaces[edge][0]] <<" , " << faceIntPoint[ edgeFaces[edge][1]]);
00073         }
00074       }
00075       // if there is more than one point than we ignore all of them, and just print a warning
00076       if (nrPoints > 1){
00077         SUNDANCE_MSG1(verb, " WARNING: nr of intersection points: " << nrPoints << " for edge [ " << brickPoints[edegIndex[edge][0]] << " , "
00078             << brickPoints[edegIndex[edge][1]] << "] , the first two points i1:" << result[0] << " , i2: " << result[1]);
00079       }
00080   }
00081 
00082   // this means no surface point has found
00083   if ( t_inters_points < 3){
00084     intersectPoints.resize(0);
00085     edgeIndex.resize(0);
00086     triangleIndex.resize(0);
00087     SUNDANCE_MSG1(verb, " WARNING: less than 3 int point found : " << t_inters_points << " for maxCellLID=" << maxCellLID
00088                 << " P[0]=" << brickPoints[0] );
00089     return;
00090   }
00091 
00092 // ------------------------ discover the triangles which form the surf -------------------------
00093   // get the polygon
00094   Array<int> polygonList(t_inters_points,-1);
00095   Array<int> pointUsed(t_inters_points,-1);
00096   int startIndex = -1 ;
00097   polygonList[0] = 0; pointUsed[0] = 1;
00098   SUNDANCE_MSG3(verb, " polygonList[ 0 ] = 0 , P: " << intersectPoints[ 0 ]);
00099   for (int jj = 1 ; jj < t_inters_points ; jj++){
00100     // look for a point which is connected to polygonList[jj] ad has not been used yet
00101     tmp_i = -1;
00102     for (int ii = 0 ; ii < t_inters_points ; ii++){
00103       // test if this edge is neighboring
00104       if ( (pointUsed[ii] < 0) && (edgeNeighboredges[ edgeIndex[polygonList[jj-1]] ][ edgeIndex[ii] ] > 0) ) {
00105         // we found one neighboring point
00106         tmp_i = ii;
00107         polygonList[jj] = ii;
00108         pointUsed[ii] = 1;
00109         SUNDANCE_MSG3(verb, " polygonList[ " << jj << " ] = " << ii << " , P: " << intersectPoints[ ii ]);
00110         break;
00111       }
00112     }
00113     TEUCHOS_TEST_FOR_EXCEPTION( tmp_i < 0 , std::runtime_error , " error by polygon detecting tmp_i > -1 , might be that the surface is not KONVEX");
00114     // look for an intersection point which is in the same face as the current intersection point
00115   }
00116 
00117   //get the starting point by getting the point which is closest
00118   Array<double> dist(t_inters_points);
00119   double dist_tmp = 1e+100;
00120   for (int jj = 0 ; jj < t_inters_points ; jj++){
00121     double sum = 0.0;
00122     for (int ii = 0 ; ii < t_inters_points ; ii++){
00123       sum = sum + ::sqrt( (intersectPoints[ polygonList[ii] ] - intersectPoints[ polygonList[jj] ])
00124           *(intersectPoints[ polygonList[ii] ] - intersectPoints[ polygonList[jj] ]) );
00125     }
00126     // store the minimum distance and the index
00127     if (sum < dist_tmp) { startIndex = jj; dist_tmp = sum; }
00128   }
00129   SUNDANCE_MSG3(verb, " central point found: " << startIndex << " P = " << intersectPoints[ polygonList[startIndex] ]);
00130 
00131 // ----------------- having the polygon and the starting point define the triangles ---------------
00132   int fI = 0 , actIF = startIndex , actIL = startIndex;
00133   triangleIndex.resize(3*(t_inters_points-2),-1);
00134   for (int jj = 0 ; jj < t_inters_points-2 ; jj++){
00135     actIF = (actIF+1) % t_inters_points;
00136     actIL = (actIF+1) % t_inters_points;
00137     SUNDANCE_MSG3(verb, " startIndex = " << startIndex << " , actIF = " << actIF << " , actIL = " << actIL);
00138     triangleIndex[fI]   = polygonList[startIndex];
00139     triangleIndex[fI+1] = polygonList[actIF];
00140     triangleIndex[fI+2] = polygonList[actIL];
00141     SUNDANCE_MSG3(verb, "add triangle [ " << polygonList[startIndex] << " , " << polygonList[actIF] << " , " << polygonList[actIL] << " ]");
00142     SUNDANCE_MSG3(verb, "add triangle points = [ " << intersectPoints[polygonList[startIndex]] << " , " << intersectPoints[polygonList[actIF]]
00143                           << " , " << intersectPoints[ polygonList[actIL] ] << " ]");
00144     fI = fI + 3;
00145   }
00146 // end the triangle is set
00147 }