00001
00002
00003
00004
00005
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
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
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
00021 const int SundanceSurf3DCalc::edgeNeighboredges[12][12] = { { 1 , 1 , 1 , 1 , 1 , 1 , 0 , 0 , 1 , 0 , 0 , 0} ,
00022 { 1 , 1 , 1 , 1 , 0 , 1 , 1 , 0 , 0 , 1 , 0 , 0} ,
00023 { 1 , 1 , 1 , 0 , 1 , 0 , 1 , 0 , 1 , 1 , 0 , 0} ,
00024 { 1 , 1 , 0 , 1 , 1 , 1 , 0 , 1 , 0 , 0 , 1 , 0} ,
00025 { 1 , 0 , 1 , 1 , 1 , 0 , 0 , 1 , 1 , 0 , 1 , 0} ,
00026 { 1 , 1 , 0 , 1 , 0 , 1 , 1 , 1 , 0 , 0 , 0 , 1} ,
00027 { 0 , 1 , 1 , 0 , 0 , 1 , 1 , 1 , 0 , 1 , 0 , 1} ,
00028 { 0 , 0 , 0 , 1 , 1 , 1 , 1 , 1 , 0 , 0 , 1 , 1} ,
00029 { 1 , 0 , 1 , 0 , 1 , 0 , 0 , 0 , 1 , 1 , 1 , 1} ,
00030 { 0 , 1 , 1 , 0 , 0 , 0 , 1 , 0 , 1 , 1 , 1 , 1} ,
00031 { 0 , 0 , 0 , 1 , 1 , 0 , 0 , 1 , 1 , 1 , 1 , 1} ,
00032 { 0 , 0 , 0 , 0 , 0 , 1 , 1 , 1 , 1 , 1 , 1 , 1} , };
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
00052
00053
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
00058 if (nrPoints == 1){
00059
00060
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
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
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
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
00093
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
00101 tmp_i = -1;
00102 for (int ii = 0 ; ii < t_inters_points ; ii++){
00103
00104 if ( (pointUsed[ii] < 0) && (edgeNeighboredges[ edgeIndex[polygonList[jj-1]] ][ edgeIndex[ii] ] > 0) ) {
00105
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 }