00001
00002
00003
00004
00005
00006
00007
00008 #include "SundanceCurveIntegralCalc.hpp"
00009 #include "SundanceOut.hpp"
00010 #include "PlayaTabs.hpp"
00011 #include "SundanceCellJacobianBatch.hpp"
00012 #include "SundancePolygon2D.hpp"
00013 #include "SundancePolygonQuadrature.hpp"
00014 #include "SundanceSurf3DCalc.hpp"
00015 #include "SundanceSurfQuadrature.hpp"
00016
00017 using namespace Sundance;
00018
00019 CurveIntegralCalc::CurveIntegralCalc() {
00020
00021 }
00022
00023 void CurveIntegralCalc::getCurveQuadPoints(CellType maxCellType ,
00024 int maxCellLID ,
00025 const Mesh& mesh ,
00026 const ParametrizedCurve& paramCurve,
00027 const QuadratureFamily& quad ,
00028 Array<Point>& curvePoints ,
00029 Array<Point>& curveDerivs ,
00030 Array<Point>& curveNormals ){
00031
00032 int verb = 0;
00033 Tabs tabs;
00034
00035
00036 int nr_point = mesh.numFacets( mesh.spatialDim() , 0 , 0);
00037 Array<Point> maxCellsPoints(nr_point);
00038 int tmp_i , point_LID;
00039
00040 SUNDANCE_MSG3(verb, tabs << "CurveIntegralCalc::getCurveQuadPoints nr points per cell: " << nr_point)
00041 for (int jj = 0 ; jj < nr_point ; jj++){
00042 point_LID = mesh.facetLID( mesh.spatialDim() , maxCellLID , 0 , jj , tmp_i );
00043 maxCellsPoints[jj] = mesh.nodePosition(point_LID);
00044 SUNDANCE_MSG3(verb, tabs << " max cell point p[" << jj << "]:"<< maxCellsPoints[jj]);
00045 }
00046
00047
00048 if (usePolygoneCurveQuadrature( maxCellType , mesh , paramCurve)){
00049
00050 CurveIntegralCalc::getCurveQuadPoints_polygon( maxCellLID , maxCellType , mesh , maxCellsPoints , paramCurve,
00051 quad , curvePoints , curveDerivs , curveNormals);
00052 } else {
00053
00054 if ( use3DSurfQuadrature( maxCellType , mesh ) ){
00055 CurveIntegralCalc::getSurfQuadPoints(maxCellLID , maxCellType , mesh , maxCellsPoints , paramCurve,
00056 quad , curvePoints , curveDerivs , curveNormals);
00057
00058
00059 } else {
00060 CurveIntegralCalc::getCurveQuadPoints_line( maxCellLID , maxCellType , mesh , maxCellsPoints , paramCurve,
00061 quad , curvePoints , curveDerivs , curveNormals);
00062 }
00063 }
00064
00065 }
00066
00067 void CurveIntegralCalc::getCurveQuadPoints_line(
00068 int maxCellLID ,
00069 CellType maxCellType ,
00070 const Mesh& mesh ,
00071 const Array<Point>& cellPoints,
00072 const ParametrizedCurve& paramCurve,
00073 const QuadratureFamily& quad ,
00074 Array<Point>& curvePoints ,
00075 Array<Point>& curveDerivs ,
00076 Array<Point>& curveNormals ){
00077
00078 int verb = 0;
00079 Tabs tabs;
00080
00081
00082 switch (paramCurve.getCurveDim()){
00083 case 1: {
00084
00085 SUNDANCE_MSG3(verb, tabs << "CurveIntegralCalc::getCurveQuadPoints_line, curve has ONE dimnesion")
00086
00087
00088
00089 Point startPoint(0.0,0.0);
00090 Point endPoint(0.0,0.0);
00091 int nrPoints , total_points = 0 ;
00092
00093 switch (maxCellType){
00094 case QuadCell:{
00095
00096
00097 TEUCHOS_TEST_FOR_EXCEPTION( cellPoints.size() != 4 ,
00098 std::runtime_error ," CurveIntegralCalc::getCurveQuadPoints , QuadCell must have 4 points" );
00099 SUNDANCE_MSG3(verb, tabs << "CurveIntegralCalc::getCurveQuadPoints_line, on QuadCell")
00100 Array<Point> result(0);
00101 int edegIndex[4][2] = { {0,1} , {0,2} , {1,3} , {2,3} };
00102
00103
00104 for (int jj = 0 ; jj < 4 ; jj++ ){
00105 paramCurve.returnIntersectPoints(cellPoints[edegIndex[jj][0]], cellPoints[edegIndex[jj][1]], nrPoints , result);
00106
00107 if (nrPoints == 1){
00108 if (total_points == 0) startPoint = result[0];
00109 else endPoint = result[0];
00110 SUNDANCE_MSG3(verb, tabs << "found Int. point:" << result[0]);
00111 total_points += nrPoints;
00112 }
00113 SUNDANCE_MSG3(verb, tabs << "ind:" << jj << ", nr Int points :" << nrPoints
00114 << " , start:" << startPoint << ", end:"<< endPoint);
00115
00116 TEUCHOS_TEST_FOR_EXCEPTION( nrPoints > 1 ,
00117 std::runtime_error , " CurveIntegralCalc::getCurveQuadPoints , QuadCell one edge " << jj << " , can have only one intersection point"
00118 << " edgeP0:" << cellPoints[edegIndex[jj][0]] << " edgeP1:" << cellPoints[edegIndex[jj][1]] );
00119 }
00120
00121 TEUCHOS_TEST_FOR_EXCEPTION( total_points > 2 ,std::runtime_error , " CurveIntegralCalc::getCurveQuadPoints total_points > 2 : " << total_points );
00122 SUNDANCE_MSG3(verb, tabs << "CurveIntegralCalc::getCurveQuadPoints_line, end finding intersection points")
00123 } break;
00124 case TriangleCell:{
00125 TEUCHOS_TEST_FOR_EXCEPTION( cellPoints.size() != 3 , std::runtime_error , " CurveIntegralCalc::getCurveQuadPoints , TriangleCell must have 3 points" );
00126 Array<Point> result;
00127 int edegIndex[3][2] = { {0,1} , {0,2} , {1,2} };
00128
00129
00130 for (int jj = 0 ; jj < 3 ; jj++ ){
00131 paramCurve.returnIntersectPoints(cellPoints[edegIndex[jj][0]], cellPoints[edegIndex[jj][1]], nrPoints , result);
00132
00133 if (nrPoints == 1){
00134 if (total_points == 0) startPoint = result[0];
00135 else endPoint = result[0];
00136 SUNDANCE_MSG3(verb, tabs << "found Int. point:" << result[0]);
00137 total_points += nrPoints;
00138 }
00139
00140 TEUCHOS_TEST_FOR_EXCEPTION( nrPoints > 1 ,
00141 std::runtime_error , " CurveIntegralCalc::getCurveQuadPoints , TriangleCell one edge " << jj << " , can have only one intersection point"
00142 << " edgeP0:" << cellPoints[edegIndex[jj][0]] << " edgeP1:" << cellPoints[edegIndex[jj][1]] );
00143 }
00144
00145 TEUCHOS_TEST_FOR_EXCEPTION( total_points > 2 ,std::runtime_error , " CurveIntegralCalc::getCurveQuadPoints total_points > 2 : " << total_points );
00146 } break;
00147 default : {
00148 TEUCHOS_TEST_FOR_EXCEPTION( true , std::runtime_error , "CurveIntegralCalc::getCurveQuadPoints , Unknown Cell in 2D" );
00149 }
00150 }
00151
00152 TEUCHOS_TEST_FOR_EXCEPTION( total_points > 2 ,std::runtime_error , " CurveIntegralCalc::getCurveQuadPoints , no much intersection points: " << total_points);
00153
00154
00155
00156
00157
00158
00159
00160 if (startPoint[0] > endPoint[0]){
00161 Point tmp = startPoint;
00162 startPoint = endPoint;
00163 endPoint = tmp;
00164 }
00165 SUNDANCE_MSG3(verb, tabs << "start end and points , start:" << startPoint << ", end:"<< endPoint)
00166
00167
00168 Array<Point> quadPoints;
00169 Array<double> quadWeights;
00170 quad.getPoints( LineCell , quadPoints , quadWeights );
00171 int nr_quad_points = quadPoints.size();
00172
00173
00174 curvePoints.resize(nr_quad_points);
00175 SUNDANCE_MSG3(verb, tabs << " setting reference quadrature points" );
00176 for (int ii = 0 ; ii < nr_quad_points ; ii++) {
00177 SUNDANCE_MSG3(verb, tabs << " nr:" << ii << " Line quad point:" << quadPoints[ii] << ", line:" << (endPoint - startPoint));
00178 curvePoints[ii] = startPoint + quadPoints[ii][0]*(endPoint - startPoint);
00179
00180
00181
00182
00183
00184
00185
00186 Array<int> cellLID(1); cellLID[0] = maxCellLID;
00187 CellJacobianBatch JBatch;
00188 JBatch.resize(1, 2, 2);
00189 mesh.getJacobians(2, cellLID , JBatch);
00190 double pointc[2] = { curvePoints[ii][0] - cellPoints[0][0] , curvePoints[ii][1] - cellPoints[0][1]};
00191 JBatch.applyInvJ(0, pointc , 1 , false);
00192 curvePoints[ii][0] = pointc[0];
00193 curvePoints[ii][1] = pointc[1];
00194
00195 SUNDANCE_MSG3(verb, tabs << " quad point nr:" << ii << " = " << curvePoints[ii]);
00196 }
00197
00198
00199
00200 curveDerivs.resize(nr_quad_points);
00201 Point dist_point(endPoint[0] - startPoint[0] , endPoint[1] - startPoint[1]);
00202 SUNDANCE_MSG3(verb, tabs << " setting derivative values points" );
00203 for (int ii = 0 ; ii < nr_quad_points ; ii++) {
00204 curveDerivs[ii] = endPoint - startPoint;
00205 SUNDANCE_MSG3(verb, tabs << " curve Derivatives point nr:" << ii << " = " << curveDerivs[ii]);
00206 }
00207
00208
00209
00210 curveNormals.resize(nr_quad_points);
00211
00212
00213 SUNDANCE_MSG3(verb, tabs << " setting normal values at points" );
00214 for (int ii = 0 ; ii < nr_quad_points ; ii++) {
00215 Point norm_vec = Point(-curveDerivs[ii][1], curveDerivs[ii][0]);
00216 norm_vec = norm_vec / ::sqrt(norm_vec*norm_vec);
00217 double relative_length = 1e-2 * ::sqrt((endPoint - startPoint)*(endPoint - startPoint));
00218 if ( paramCurve.curveEquation( startPoint + relative_length*norm_vec ) > 0 ) {
00219 curveNormals[ii] = -norm_vec;
00220 }
00221 else {
00222 curveNormals[ii] = norm_vec;
00223 }
00224 SUNDANCE_MSG3(verb, tabs << " curve Normals point nr:" << ii << " = " << curveNormals[ii]);
00225 }
00226
00227
00228 } break;
00229
00230
00231
00232
00233
00234 case 2: {
00235
00236
00237
00238
00239
00240
00241
00242
00243 Tabs tabs;
00244 switch (maxCellType){
00245 case BrickCell:{
00246 int edegIndex[12][2] = { {0,1} , {0,2} , {0,4} , {1,3} , {1,5} , {2,3} , {2,6} , {3,7} ,
00247 {4,5} , {4,6} , {5,7} , {6,7} };
00248 int 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} };
00249
00250 int t_inters_points = 0;
00251 int nrPoints = 0;
00252 Array<Point> edgeIntersectPoint(t_inters_points);
00253 Array<int> edgeIndex(t_inters_points);
00254 Array<int> edgePointI(12,-1);
00255
00256
00257 for (int jj = 0 ; jj < 12 ; jj++ ){
00258 Array<Point> result(0);
00259 paramCurve.returnIntersectPoints(cellPoints[edegIndex[jj][0]], cellPoints[edegIndex[jj][1]], nrPoints , result);
00260
00261 if (nrPoints > 0){
00262 edgeIntersectPoint.resize( t_inters_points + 1 );
00263 edgeIndex.resize( t_inters_points + 1 );
00264 edgeIntersectPoint[t_inters_points] = result[0];
00265 edgeIndex[t_inters_points] = jj;
00266 edgePointI[jj] = t_inters_points;
00267 SUNDANCE_MSG3(verb, tabs << "found Int. point [" << t_inters_points <<"]=" << result[0]);
00268 t_inters_points++;
00269 }
00270 SUNDANCE_MSG3(verb, tabs << "ind:" << jj << ", nr Int points :" << nrPoints );
00271 TEUCHOS_TEST_FOR_EXCEPTION( nrPoints > 1 ,
00272 std::runtime_error , " CurveIntegralCalc::getCurveQuadPoints , QuadCell one edge " << jj << " , can have only one intersection point" );
00273 }
00274
00275
00276
00277 Array<double> totDist(t_inters_points);
00278
00279 SUNDANCE_MSG3(verb, tabs << " Calculating distances ");
00280 for (int i = 0 ; i < t_inters_points ; i++){
00281 double sum = 0.0;
00282 for (int j = 0 ; j < t_inters_points ; j++){
00283 sum = sum + ::sqrt(edgeIntersectPoint[i]*edgeIntersectPoint[j]);
00284 }
00285 totDist[i] = sum;
00286 }
00287
00288 int minIndex = 0;
00289 for (int i = 1 ; i < t_inters_points ; i++){
00290 if (totDist[i] < totDist[minIndex]) minIndex = i;
00291 }
00292 SUNDANCE_MSG3(verb, tabs << " minIndex = " << minIndex );
00293
00294
00295
00296
00297 int intersTriangleInd = 0;
00298 int nrTriangles = t_inters_points-2;
00299 Array<int> triangles(3*nrTriangles,-1);
00300 Array<double> triangle_area(nrTriangles,0.0);
00301 double total_area = 0.0;
00302 SUNDANCE_MSG3(verb, tabs << " nrTriangles = " << nrTriangles );
00303 int intPointProc = -1;
00304 int intPointProc_lastNeigh = -1;
00305
00306
00307
00308 for (int jj = 0 ; jj < 6 ; jj++ )
00309 {
00310 int nrIpointPerFace = 0;
00311 int firstPI = -1 , secondPI = -1;
00312
00313 for (int ii = 0 ; ii < 4 ; ii++){
00314
00315 if ( edgePointI[faceEdges[jj][ii]] >= 0){
00316 if (nrIpointPerFace > 0){
00317 secondPI = edgePointI[faceEdges[jj][ii]];
00318 }else{
00319 firstPI = edgePointI[faceEdges[jj][ii]];
00320 }
00321 TEUCHOS_TEST_FOR_EXCEPTION( nrIpointPerFace > 2 , std::runtime_error,"getCurveQuadPoints , nrIpointPerFace > 2 , " << nrIpointPerFace );
00322 nrIpointPerFace++;
00323 }
00324 }
00325 TEUCHOS_TEST_FOR_EXCEPTION( nrIpointPerFace == 1 , std::runtime_error,"getCurveQuadPoints , nrIpointPerFace == 1 , " << nrIpointPerFace );
00326
00327 if (( nrIpointPerFace > 1) && ((minIndex == firstPI) || (minIndex == secondPI) )){
00328 SUNDANCE_MSG3(verb, tabs << " Found starting line:" << firstPI << " -> " << secondPI);
00329 intPointProc_lastNeigh = minIndex;
00330 if (minIndex == firstPI){
00331 intPointProc = secondPI;
00332 }else{
00333 intPointProc = firstPI;
00334 }
00335
00336 break;
00337 }
00338 }
00339 TEUCHOS_TEST_FOR_EXCEPTION( intPointProc < 0 , std::runtime_error,"getCurveQuadPoints , intPointProc < 0 , " << intPointProc );
00340
00341
00342 for (int pI = 0 ; pI < nrTriangles; pI++)
00343 {
00344
00345 for (int jj = 0 ; jj < 6 ; jj++ )
00346 {
00347 int nrIpointPerFace = 0;
00348 int firstPI = -1 , secondPI = -1;
00349
00350 for (int ii = 0 ; ii < 4 ; ii++){
00351
00352 if ( edgePointI[faceEdges[jj][ii]] >= 0){
00353 if (nrIpointPerFace > 0){
00354 secondPI = edgePointI[faceEdges[jj][ii]];
00355 }else{
00356 firstPI = edgePointI[faceEdges[jj][ii]];
00357 }
00358 TEUCHOS_TEST_FOR_EXCEPTION( nrIpointPerFace > 2 , std::runtime_error,"getCurveQuadPoints , nrIpointPerFace > 2 , " << nrIpointPerFace );
00359 nrIpointPerFace++;
00360 }
00361 }
00362 TEUCHOS_TEST_FOR_EXCEPTION( nrIpointPerFace == 1 , std::runtime_error,"getCurveQuadPoints , nrIpointPerFace == 1 , " << nrIpointPerFace );
00363
00364 if (( nrIpointPerFace > 1) && ((intPointProc == firstPI) || (intPointProc == secondPI))
00365 && ((intPointProc_lastNeigh != firstPI) && (intPointProc_lastNeigh != secondPI)))
00366 {
00367 SUNDANCE_MSG3(verb, tabs << " Found next line:" << firstPI << " -> " << secondPI);
00368 if (intPointProc == firstPI){
00369 intPointProc = secondPI;
00370 intPointProc_lastNeigh = firstPI;
00371 }else{
00372 intPointProc = firstPI;
00373 intPointProc_lastNeigh = secondPI;
00374 }
00375
00376 triangles[intersTriangleInd] = minIndex;
00377 triangles[intersTriangleInd+1] = intPointProc_lastNeigh;
00378 triangles[intersTriangleInd+2] = intPointProc;
00379 SUNDANCE_MSG3(verb, tabs << " Found triangle:" << minIndex << " , " << firstPI << " , " << secondPI);
00380 Point v1 = edgeIntersectPoint[firstPI] - edgeIntersectPoint[minIndex];
00381 Point v2 = edgeIntersectPoint[secondPI] - edgeIntersectPoint[minIndex];
00382 Point areaV = cross(v1,v2);
00383 SUNDANCE_MSG3(verb, tabs << " TriangINdex = " << intersTriangleInd << " , area = " << ::sqrt(areaV*areaV));
00384 triangle_area[ intersTriangleInd/3 ] = ::sqrt(areaV*areaV) / (double)2.0 ;
00385 total_area = total_area + triangle_area[ intersTriangleInd/3 ];
00386 intersTriangleInd = intersTriangleInd + 3;
00387
00388 break;
00389 }
00390 }
00391 }
00392
00393 SUNDANCE_MSG3(verb, tabs << " END Found triangle total_area=" << total_area);
00394
00395
00396
00397
00398 if (t_inters_points < 4){
00399
00400 SUNDANCE_MSG3(verb, tabs << " Add one additional point to have at least 2 triangles ");
00401 int longEdI1 = -1 , longEdI2 = -1;
00402 for (int ii = 0 ; ii < t_inters_points ; ii++)
00403 if ( (ii != minIndex) ){
00404 longEdI1 = ii;
00405 break;
00406 }
00407 for (int ii = 0 ; ii < t_inters_points ; ii++)
00408 if ( (ii != minIndex) && (ii != longEdI1) ){
00409 longEdI2 = ii;
00410 break;
00411 }
00412 TEUCHOS_TEST_FOR_EXCEPTION( (longEdI1 < 0)||(longEdI2 < 0), std::runtime_error," (longEdI1 < 0)||(longEdI2 < 0):" << longEdI1 << "," <<longEdI2);
00413
00414
00415 edgeIntersectPoint.resize(t_inters_points+1);
00416 edgeIntersectPoint[t_inters_points] = edgeIntersectPoint[longEdI1]/2.0 + edgeIntersectPoint[longEdI2]/2.0;
00417
00418 int new_p = t_inters_points;
00419 t_inters_points += 1;
00420 nrTriangles += 1;
00421
00422
00423 triangles.resize(3*(nrTriangles));
00424 triangles[0] = minIndex;
00425 triangles[1] = longEdI1;
00426 triangles[2] = new_p;
00427 triangles[3] = minIndex;
00428 triangles[4] = new_p;
00429 triangles[5] = longEdI2;
00430
00431 triangle_area[0] = triangle_area[0] / 2.0;
00432 triangle_area.resize(2);
00433 triangle_area[1] = triangle_area[0];
00434 }
00435
00436
00437 Array<Point> quadPoints;
00438 Array<double> quadWeights;
00439 quad.getPoints( QuadCell , quadPoints , quadWeights );
00440 int nr_quad_points = quadPoints.size();
00441
00442
00443 curvePoints.resize(nr_quad_points);
00444 Array<int> triangleInd(nr_quad_points,-1);
00445 SUNDANCE_MSG3(verb, tabs << " setting reference quadrature points" );
00446 for (int ii = 0 ; ii < nr_quad_points ; ii++) {
00447
00448 int tInd = -1;
00449 Point tmp(0.0,0.0,0.0);
00450 get3DRealCoordsOnSurf( quadPoints[ii] , cellPoints , triangles ,
00451 nrTriangles , edgeIntersectPoint, tInd , tmp );
00452 triangleInd[ii] = tInd;
00453 curvePoints[ii] = tmp;
00454
00455 SUNDANCE_MSG3(verb, tabs << " quad point nr:" << ii << " = " << curvePoints[ii]);
00456 }
00457
00458
00459 curveDerivs.resize(nr_quad_points);
00460 SUNDANCE_MSG3(verb, tabs << " setting surface values points" );
00461 for (int ii = 0 ; ii < nr_quad_points ; ii++) {
00462
00463 int tInd = triangleInd[ii] ;
00464 curveDerivs[ii] = Point( nrTriangles*triangle_area[tInd] , 0.0, 0.0 );
00465
00466 SUNDANCE_MSG3(verb, tabs << " curve Derivatives point nr:" << ii << " = " << curveDerivs[ii]);
00467 }
00468
00469
00470
00471 curveNormals.resize(nr_quad_points);
00472 SUNDANCE_MSG3(verb, tabs << " setting normal values at points" );
00473 for (int ii = 0 ; ii < nr_quad_points ; ii++) {
00474 int tInd = triangleInd[ii];
00475
00476 Point norm_vec = cross( edgeIntersectPoint[triangles[3*tInd+1]] - edgeIntersectPoint[triangles[3*tInd]] ,
00477 edgeIntersectPoint[triangles[3*tInd+2]] - edgeIntersectPoint[triangles[3*tInd]]);
00478 double relative_length = 1e-2 * ::sqrt(norm_vec*norm_vec);
00479 norm_vec = norm_vec / ::sqrt(norm_vec*norm_vec);
00480 if ( paramCurve.curveEquation( edgeIntersectPoint[triangles[3*tInd]] + relative_length*norm_vec ) > 0 ) {
00481 curveNormals[ii] = -norm_vec;
00482 }
00483 else {
00484 curveNormals[ii] = norm_vec;
00485 }
00486 SUNDANCE_MSG3(verb, tabs << " curve Normals point nr:" << ii << " = " << curveNormals[ii]);
00487 }
00488
00489 break;
00490 }
00491 case TetCell:{
00492
00493 TEUCHOS_TEST_FOR_EXCEPTION( true, std::runtime_error,"CurveIntagralCalc::getCurveQuadPoints , not implemented for " << maxCellType );
00494 break;
00495 }
00496 default:{
00497 TEUCHOS_TEST_FOR_EXCEPTION( true, std::runtime_error,"CurveIntagralCalc::getCurveQuadPoints , not implemented for " << maxCellType );
00498 }
00499 }
00500 SUNDANCE_MSG3(verb, tabs << " END CurveIntagralCalc::getCurveQuadPoints");
00501 } break;
00502
00503
00504
00505 default: {
00506
00507 TEUCHOS_TEST_FOR_EXCEPTION( true, std::runtime_error,"CurveIntagralCalc::getCurveQuadPoints , curve dimension must be 1 or two ");
00508 }
00509 }
00510 }
00511
00512
00513 void CurveIntegralCalc::get3DRealCoordsOnSurf(const Point &refP ,
00514 const Array<Point>& cellPoints,
00515 const Array<int> &triangles ,
00516 const int nrTriag ,
00517 const Array<Point> edgeIntersectPoint,
00518 int &triagIndex ,
00519 Point &realPoint){
00520
00521 Tabs tabs;
00522 int verb = 0;
00523 realPoint[0] = 0.0; realPoint[1] = 0.0; realPoint[2] = 0.0;
00524 SUNDANCE_MSG3(verb, tabs << " start get3DRealCoordsOnSurf refP="<<refP );
00525 Point v1;
00526 Point v2;
00527
00528 double case1[2][4] = {{1 ,-1 ,0 ,1},{ 1 , 0 , -1 , 1}};
00529 double case2[3][4] = {{1 ,-2, 0, 2} , { 2 , -2, -1, 2},{ 1 , 0, -1, 1}};
00530 double case3[4][4] = {{1, -2,0 , 2} , { 2 , -2, -1,2},{ 2, -1,-2, 2},{ 2 ,0, -2, 1}};
00531 double *refInv = 0;
00532 switch (nrTriag){
00533 case 2:{
00534 if (refP[0] >= refP[1]){
00535 triagIndex = 0; v1 = Point(1.0,0.0); v2 = Point(1.0,1.0);
00536 }else{
00537 triagIndex = 1; v1 = Point(1.0,1.0); v2 = Point(0.0,1.0);
00538 }
00539 refInv = case1[triagIndex];
00540 break;
00541 }
00542 case 3:{
00543 if (refP[0] >= 2*refP[1]){
00544 triagIndex = 0; v1 = Point(1.0,0.0); v2 = Point(1.0,0.5);
00545 }else if ((refP[0] <= 2*refP[1]) && (refP[0] >= refP[1])){
00546 triagIndex = 1; v1 = Point(1.0,0.5); v2 = Point(1.0,1.0);
00547 }else{
00548 triagIndex = 2; v1 = Point(1.0,1.0); v2 = Point(0.0,1.0);
00549 }
00550 refInv = case2[triagIndex];
00551 break;
00552 }
00553 case 4:{
00554 if (refP[0] >= 2*refP[1]){
00555 triagIndex = 0; v1 = Point(1.0,0.0); v2 = Point(1.0,0.5);
00556 }else if ((refP[0] <= 2*refP[1]) && (refP[0] >= refP[1])){
00557 triagIndex = 1; v1 = Point(1.0,0.5); v2 = Point(1.0,1.0);
00558 }else if ((refP[0] <= refP[1]) && (2*refP[0] >= refP[1])){
00559 triagIndex = 2; v1 = Point(1.0,1.0); v2 = Point(0.5,1.0);
00560 }else{
00561 triagIndex = 3;v1 = Point(0.5,1.0); v2 = Point(0.0,1.0);
00562 }
00563 refInv = case3[triagIndex];
00564 break;
00565 }
00566 default:{
00567 TEUCHOS_TEST_FOR_EXCEPTION( true, std::runtime_error,"get3DRealCoordsOnSurf , too many triangles " << nrTriag);
00568 }
00569 }
00570
00571 Point p0 = edgeIntersectPoint[triangles[3*triagIndex]];
00572 Point p1 = edgeIntersectPoint[triangles[3*triagIndex+1]];
00573 Point p2 = edgeIntersectPoint[triangles[3*triagIndex+2]];
00574 SUNDANCE_MSG3(verb, tabs << "get3DRealCoordsOnSurf p0="<<p0<<" , p1="<<p1<<" , p2="<<p2);
00575 SUNDANCE_MSG3(verb, tabs << "get3DRealCoordsOnSurf v1=" << v1 << " , v2=" << v2 << " , triagIndex=" << triagIndex);
00576
00577
00578 double ref_x = refInv[0] * refP[0] + refInv[1]*refP[1];
00579 double ref_y = refInv[2] * refP[0] + refInv[3]*refP[1];
00580 SUNDANCE_MSG3(verb, tabs << " after traf to ref_x ="<<ref_x<<" , ref_y="<<ref_y );
00581
00582
00583 realPoint = p0 + ref_x*(p1-p0) + ref_y*(p2-p0);
00584 SUNDANCE_MSG3(verb, tabs << " end get3DRealCoordsOnSurf , realPoint="<<realPoint );
00585 SUNDANCE_MSG3(verb, tabs << " end get3DRealCoordsOnSurf cellPoints[0]="<<cellPoints[0]<<" , cellPoints[7]="<<cellPoints[7] );
00586
00587
00588
00589 realPoint[0] = (realPoint[0] - cellPoints[0][0]) / (cellPoints[7][0] - cellPoints[0][0]);
00590 realPoint[1] = (realPoint[1] - cellPoints[0][1]) / (cellPoints[7][1] - cellPoints[0][1]);
00591 realPoint[2] = (realPoint[2] - cellPoints[0][2]) / (cellPoints[7][2] - cellPoints[0][2]);
00592 SUNDANCE_MSG3(verb, tabs << " end get3DRealCoordsOnSurf triagIndex="<<triagIndex<<" , realPoint="<<realPoint );
00593
00594 }
00595
00596
00597 bool CurveIntegralCalc::usePolygoneCurveQuadrature(
00598 const CellType& cellType ,
00599 const Mesh& mesh ,
00600 const ParametrizedCurve& paramCurve){
00601
00602
00603 CellType maxCellType = mesh.cellType(mesh.spatialDim());
00604
00605 switch (maxCellType)
00606 {
00607
00608 case TriangleCell:
00609 {
00610 switch (cellType)
00611 {
00612 case TriangleCell:
00613 { break; }
00614 default: {}
00615 }
00616 break;
00617 }
00618 case QuadCell:
00619 {
00620 switch(cellType)
00621 {
00622 case QuadCell:
00623 {
00624
00625 const Polygon2D* poly = dynamic_cast<const Polygon2D* > ( paramCurve.ptr().get()->getUnderlyingCurveObj() );
00626 if ( poly != 0 ){
00627 return true;
00628 }
00629 break;
00630 }
00631 default: {}
00632 }
00633 break;
00634 }
00635 default: { }
00636 }
00637
00638 return false;
00639 }
00640
00641
00642 void CurveIntegralCalc::getCurveQuadPoints_polygon(
00643 int maxCellLID ,
00644 CellType maxCellType ,
00645 const Mesh& mesh ,
00646 const Array<Point>& cellPoints,
00647 const ParametrizedCurve& paramCurve,
00648 const QuadratureFamily& quad ,
00649 Array<Point>& curvePoints ,
00650 Array<Point>& curveDerivs ,
00651 Array<Point>& curveNormals ) {
00652
00653
00654 const Polygon2D* polygon = dynamic_cast<const Polygon2D* > ( paramCurve.ptr().get()->getUnderlyingCurveObj() );
00655 if ( (polygon == 0) || (maxCellType != QuadCell) ) {
00656 TEUCHOS_TEST_FOR_EXCEPTION( true , std::runtime_error , " getCurveQuadPoints_polygon , paramCurve must be a Polygon and the cell must be a quad" );
00657 return;
00658 }
00659
00660
00661
00662
00663
00664 Array<Point> allIntPoints(0);
00665 Array<int> allIntPointsEdge(0);
00666 Array<Point> allPolygonPoints(0);
00667 int nrPoints , total_points , polyNrPoint ;
00668 Tabs tabs;
00669 int verb = 0;
00670 double eps = 1e-14;
00671
00672 SUNDANCE_MSG3(verb, " ---------- START METHOD ----------- ");
00673
00674
00675 switch (maxCellType){
00676 case QuadCell:{
00677
00678
00679 TEUCHOS_TEST_FOR_EXCEPTION( cellPoints.size() != 4 ,
00680 std::runtime_error ," CurveIntegralCalc::getCurveQuadPoints , QuadCell must have 4 points" );
00681
00682 Array<Point> result(5);
00683 int edegIndex[4][2] = { {0,1} , {0,2} , {1,3} , {2,3} };
00684 bool pointExists;
00685
00686 total_points = 0;
00687
00688 for (int jj = 0 ; jj < 4 ; jj++ ){
00689 paramCurve.returnIntersectPoints(cellPoints[edegIndex[jj][0]], cellPoints[edegIndex[jj][1]], nrPoints , result);
00690
00691 for (int ptmp = 0 ; ptmp < nrPoints ; ptmp++){
00692 pointExists = false;
00693
00694 for (int tmp_r = 0 ; tmp_r < total_points ; tmp_r++){
00695 Point tmpPoint = result[ptmp];
00696 tmpPoint = tmpPoint - allIntPoints[tmp_r];
00697 if ( ::sqrt(tmpPoint * tmpPoint) < eps ) pointExists = true;
00698 }
00699
00700 if (!pointExists){
00701
00702 allIntPoints.resize( total_points + 1 );
00703 allIntPointsEdge.resize( total_points + 1 );
00704 allIntPoints[total_points] = result[ptmp];
00705 allIntPointsEdge[total_points] = jj;
00706 total_points = total_points + 1;
00707 }
00708 }
00709 }
00710
00711
00712 } break;
00713 case TriangleCell:{
00714 TEUCHOS_TEST_FOR_EXCEPTION( true , std::runtime_error , " CurveIntegralCalc::getCurveQuadPoints_polygon , TriangleCell not implemented yet" );
00715 } break;
00716 default : {
00717 TEUCHOS_TEST_FOR_EXCEPTION( true , std::runtime_error , "CurveIntegralCalc::getCurveQuadPoints_polygon , Unknown Cell in 2D" );
00718 }
00719 }
00720
00721
00722
00723 polygon->getCellsPolygonesPoint( maxCellLID ,allPolygonPoints);
00724 polyNrPoint = allPolygonPoints.size();
00725
00726
00727
00728
00729 nrPoints = polyNrPoint + total_points;
00730 Array<Point> segmentStart(0);
00731 Array<Point> segmentEnd(0);
00732 int segmentNr = 0 , flagMode , edgeActual = 0;
00733 Array<bool> usedIntPoint( total_points , false );
00734
00735
00736 if (total_points < 2) {
00737 std::cout << "total_points = " << total_points << " , P0:" << allIntPoints[0] << std::endl;
00738 }
00739 if (total_points == 3) {
00740 std::cout << "total_points = " << total_points << "P0:" << allIntPoints[0] <<
00741 " , P1:" << allIntPoints[1] << " , P2:" << allIntPoints[2] << std::endl;
00742 total_points = 2;
00743 }
00744 if (total_points > 4) {
00745 std::cout << "total_points = " << total_points << "P0:" << allIntPoints[0] << " , P1:" << allIntPoints[1]
00746 << " , P2:" << allIntPoints[2] << " , P3:" << allIntPoints[3] << std::endl;
00747 total_points = 4;
00748 }
00749
00750 TEUCHOS_TEST_FOR_EXCEPTION( (total_points != 2) && (total_points != 4) , std::runtime_error , "nr Intersect point can be eighter 2 or 4 but is " << total_points);
00751
00752
00753 flagMode = 0;
00754 for (int intP = 0 ; intP < total_points ; intP++ ){
00755 if (flagMode == 0){
00756
00757 int Pind = -1;
00758 for (int stmp = 0 ; stmp < total_points ; stmp++){
00759
00760 if (!usedIntPoint[stmp]) { Pind = stmp; break;}
00761 }
00762
00763 segmentStart.resize( segmentNr + 1);
00764 segmentEnd.resize( segmentNr + 1 );
00765 edgeActual = allIntPointsEdge[Pind];
00766 segmentStart[ segmentNr ] = allIntPoints[Pind];
00767 usedIntPoint[ Pind ] = true;
00768 flagMode = 1;
00769 } else {
00770
00771 double dist = 1e+100;
00772 int Pind = -1;
00773 for (int stmp = 0 ; stmp < total_points ; stmp++){
00774
00775 if ( !usedIntPoint[stmp] && ( (edgeActual != allIntPointsEdge[stmp]) || total_points < 3 ) ){
00776 Point tmpPointC = (segmentStart[segmentNr] - allIntPoints[stmp]);
00777 double thisDist = ::sqrt(tmpPointC*tmpPointC);
00778 if ( thisDist < dist ) {
00779 dist = thisDist;
00780 Pind = stmp;
00781 }
00782 }
00783 }
00784
00785 segmentEnd[ segmentNr ] = allIntPoints[Pind];
00786 usedIntPoint[ Pind ] = true;
00787 flagMode = 0;
00788 segmentNr = segmentNr + 1;
00789 }
00790 }
00791
00792
00793
00794 for (int polyPoint = 0 ; polyPoint < polyNrPoint ; polyPoint++ ){
00795
00796 double totalDist = 1e+100;
00797 int segInd = -1;
00798 for (int pl = 0; pl < segmentNr ; pl++ ){
00799 Point d1p = segmentStart[pl] - allPolygonPoints[polyPoint];
00800 double d1 = ::sqrt(d1p*d1p);
00801 Point d2p = segmentEnd[pl] - allPolygonPoints[polyPoint];
00802 double d2 = ::sqrt(d2p*d2p);
00803 Point dist_prev = segmentEnd[pl] - segmentStart[pl];
00804 double d_p = ::sqrt(dist_prev*dist_prev);
00805
00806
00807
00808 if ( (d1 + d2 - d_p < totalDist) && (d1 > eps) && (d2 > eps) ) {
00809 totalDist = d1 + d2 - d_p;
00810 segInd = pl;
00811 }
00812 }
00813
00814
00815
00816
00817 if ( segInd < 0) continue;
00818
00819
00820 segmentStart.resize( segmentNr + 1);
00821 segmentEnd.resize( segmentNr + 1 );
00822
00823 for (int tmpP = segmentNr-1; tmpP >= segInd ; tmpP--){
00824 segmentStart[tmpP+1] = segmentStart[tmpP];
00825 segmentEnd[tmpP+1] = segmentEnd[tmpP];
00826 }
00827
00828
00829 segmentEnd[segInd] = allPolygonPoints[polyPoint];
00830 segmentStart[segInd+1] = allPolygonPoints[polyPoint];
00831 segmentNr = segmentNr + 1;
00832 }
00833
00834
00835 Array<double> segmentLengthRel(0);
00836 Array<double> segmentStartLengthRel(0);
00837 segmentLengthRel.resize(segmentNr);
00838 segmentStartLengthRel.resize(segmentNr);
00839
00840
00841 double totalLength = 0.0 , actLen = 0.0;
00842 for (int seg = 0 ; seg < segmentNr ; seg++){
00843 totalLength = totalLength + ::sqrt( (segmentEnd[seg]-segmentStart[seg])*(segmentEnd[seg]-segmentStart[seg]) );
00844
00845
00846 }
00847
00848
00849 for (int seg = 0 ; seg < segmentNr ; seg++){
00850 double thisLength = ::sqrt( (segmentEnd[seg]-segmentStart[seg])*(segmentEnd[seg]-segmentStart[seg]) );
00851 segmentLengthRel[seg] = thisLength/totalLength;
00852 segmentStartLengthRel[seg] = actLen/totalLength;
00853
00854
00855
00856 actLen = actLen + thisLength;
00857 }
00858
00859
00860
00861 const PolygonQuadrature* polyquad = dynamic_cast<const PolygonQuadrature*> (quad.ptr().get());
00862 if ( polyquad == 0 ) {
00863 TEUCHOS_TEST_FOR_EXCEPTION( true , std::runtime_error ,
00864 " getCurveQuadPoints_polygon , Quadrature for a polygon curve must be PolygonQuadrature" );
00865 return;
00866 }
00867
00868
00869 Array<Point> quadPoints;
00870 Array<double> quadWeights;
00871 polyquad->getPoints( LineCell , quadPoints , quadWeights );
00872 int nr_quad_points = quadPoints.size();
00873 Point endPoint , startPoint;
00874 int maxLineSegments = PolygonQuadrature::getNrMaxLinePerCell();
00875 Array<int> cellLID(1); cellLID[0] = maxCellLID;
00876 CellJacobianBatch JBatch;
00877 JBatch.resize(1, 2, 2);
00878 mesh.getJacobians(2, cellLID , JBatch);
00879
00880 TEUCHOS_TEST_FOR_EXCEPTION( maxLineSegments < segmentNr , std::runtime_error , " getCurveQuadPoints_polygon , nr of polygon lines inside one cell too high" <<
00881 " Use PolygonQuadrature::getNrMaxLinePerCell(" << segmentNr << ") , now it is only: " << maxLineSegments );
00882
00883
00884
00885 curvePoints.resize(nr_quad_points);
00886 curveDerivs.resize(nr_quad_points);
00887 curveNormals.resize(nr_quad_points);
00888
00889
00890 int ii = 0;
00891 for (int seg = 0; seg < segmentNr ; seg++ )
00892 {
00893
00894 startPoint = segmentStart[seg];
00895 endPoint = segmentEnd[seg];
00896
00897
00898
00899 Point curvderi = endPoint - startPoint;
00900 Point norm_vec = Point( -curvderi[1] , curvderi[0] );
00901 norm_vec = norm_vec / ::sqrt(norm_vec*norm_vec);
00902 double relative_length = 1e-2 * ::sqrt(curvderi*curvderi);
00903 bool invNormVect = false;
00904
00905 if ( paramCurve.curveEquation( (0.5*endPoint + 0.5*startPoint) + relative_length*norm_vec ) > 0 ) {
00906 invNormVect = true;
00907 }
00908 else {
00909 invNormVect = false;
00910 }
00911
00912
00913 for (int pointNr = 0 ; pointNr < nr_quad_points/maxLineSegments ; pointNr++ , ii++ )
00914 {
00915
00916
00917
00918
00919
00920 curvePoints[ii] = startPoint + quadPoints[ii][0]*(endPoint - startPoint);
00921
00922
00923 double pointc[2] = { curvePoints[ii][0] - cellPoints[0][0] , curvePoints[ii][1] - cellPoints[0][1] };
00924 JBatch.applyInvJ(0, pointc , 1 , false);
00925 curvePoints[ii][0] = pointc[0];
00926 curvePoints[ii][1] = pointc[1];
00927
00928
00929
00930
00931 curveDerivs[ii] = Point( segmentLengthRel[seg]*totalLength*((double)maxLineSegments) , 0.0 );
00932
00933
00934
00935 if ( invNormVect ) {
00936 curveNormals[ii] = -norm_vec;
00937 }
00938 else {
00939 curveNormals[ii] = norm_vec;
00940 }
00941
00942 }
00943 }
00944
00945
00946 for ( ; ii < nr_quad_points ; ii++ ){
00947 curvePoints[ii] = Point( 0.0 , 0.0 );
00948 curveDerivs[ii] = Point( 0.0 , 0.0 );
00949 curveNormals[ii] = Point( 0.0 , 0.0 );
00950 }
00951
00952 SUNDANCE_MSG3(verb, " ---------- END METHOD ----------- ");
00953 }
00954
00955
00956 void CurveIntegralCalc::getSurfQuadPoints(
00957 int maxCellLID ,
00958 CellType maxCellType ,
00959 const Mesh& mesh ,
00960 const Array<Point>& cellPoints,
00961 const ParametrizedCurve& paramCurve,
00962 const QuadratureFamily& quad ,
00963 Array<Point>& curvePoints ,
00964 Array<Point>& curveDerivs ,
00965 Array<Point>& curveNormals ) {
00966
00967 const SurfQuadrature* surfquad = dynamic_cast<const SurfQuadrature*> (quad.ptr().get());
00968 if ( surfquad == 0 ) {
00969 TEUCHOS_TEST_FOR_EXCEPTION( true , std::runtime_error ,
00970 " getSurfQuadPoints , Surface Quadrature for a 3D must be SurfQuadrature" );
00971 return;
00972 }
00973
00974 Array<Point> quadPoints;
00975 Array<double> quadWeights;
00976 surfquad->getPoints( TriangleCell , quadPoints , quadWeights );
00977 int nr_quad_points = quadPoints.size();
00978 int maxTriangles = SurfQuadrature::getNrMaxTrianglePerCell();
00979 int nrQuadPointPerTriag = nr_quad_points/maxTriangles;
00980 int verb = 0;
00981 curvePoints.resize(nr_quad_points);
00982 curveDerivs.resize(nr_quad_points);
00983 curveNormals.resize(nr_quad_points);
00984
00985 Array<Point> intersectPoints;
00986 Array<int> edgeIndex;
00987 Array<int> triangleIndex;
00988
00989 SundanceSurf3DCalc::getCurveQuadPoints( maxCellType , maxCellLID , mesh , paramCurve, cellPoints,
00990 intersectPoints , edgeIndex , triangleIndex );
00991
00992 int nr_triag = triangleIndex.size()/3 , qind = 0;
00993 double totalArea = 0.0 , area;
00994
00995 for (int t = 0 ; t < nr_triag ; t++){
00996
00997 Point p0 = intersectPoints[triangleIndex[3*t]];
00998 Point p1 = intersectPoints[triangleIndex[3*t + 1]];
00999 Point p2 = intersectPoints[triangleIndex[3*t + 2]];
01000 Point n = cross(p1-p0,p2-p0);
01001 totalArea = totalArea + ::sqrt(n*n);
01002 }
01003
01004
01005 SUNDANCE_MSG3(verb, " nr_triag = " << nr_triag);
01006
01007 for (int t = 0 ; t < nr_triag ; t++){
01008
01009 Point p0 = intersectPoints[triangleIndex[3*t]];
01010 Point p1 = intersectPoints[triangleIndex[3*t + 1]];
01011 Point p2 = intersectPoints[triangleIndex[3*t + 2]];
01012
01013 Point n = cross(p1-p0,p2-p0);
01014
01015 area = ::sqrt(n*n);
01016
01017 if ( paramCurve.curveEquation( p0 + 5e-3*n ) > 0 ) { n = (-1)*n; }
01018 n = (1/::sqrt(n*n))*n;
01019
01020
01021
01022
01023 for (int q = 0 ; q < nrQuadPointPerTriag ; q++ , qind++){
01024
01025 Point pTmp = p0 + quadPoints[qind][0] * (p1-p0) + quadPoints[qind][1] * (p2-p0);
01026
01027
01028 curvePoints[qind] = Point( (pTmp[0] - cellPoints[0][0])/(cellPoints[7][0] - cellPoints[0][0]),
01029 (pTmp[1] - cellPoints[0][1])/(cellPoints[7][1] - cellPoints[0][1]),
01030 (pTmp[2] - cellPoints[0][2])/(cellPoints[7][2] - cellPoints[0][2]) );
01031
01032 curveDerivs[qind] = Point( area * ((double)maxTriangles) , 0.0 , 0.0 );
01033 curveNormals[qind] = n;
01034 }
01035 }
01036
01037
01038 for ( ; qind < nr_quad_points ; qind++){
01039 curvePoints[qind] = Point( 0 , 0 , 0 );
01040 curveDerivs[qind] = Point( 0 , 0 , 0 );
01041 curveNormals[qind] = Point( 0 , 0 , 0 );
01042 }
01043 SUNDANCE_MSG3(verb, " ---------- END METHOD ----------- ");
01044 }