src/meshTools/searchableSurface/searchableSurfacesQueries.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software; you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM; if not, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*----------------------------------------------------------------------------*/
  26
  27 #include "searchableSurfacesQueries.H"
  28 #include "ListOps.H"
  29 #include "OFstream.H"
  30 #include "meshTools.H"
  31
  32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  33
  34 defineTypeNameAndDebug(Foam::searchableSurfacesQueries, 0);
  35
  36 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  37
  38 Foam::pointIndexHit Foam::searchableSurfacesQueries::tempFindNearest
  39 (
  40     const searchableSurface& surf,
  41     const point& pt,
  42     const scalar initDistSqr
  43 )
  44 {
  45     pointField onePoint(1, pt);
  46     scalarField oneDist(1, initDistSqr);
  47     List<pointIndexHit> oneHit(1);
  48     surf.findNearest(onePoint, oneDist, oneHit);
  49     return oneHit[0];
  50 }
  51
  52
  53 // Calculate sum of distance to surfaces.
  54 Foam::scalar Foam::searchableSurfacesQueries::sumDistSqr
  55 (
  56     const PtrList<searchableSurface>& allSurfaces,
  57     const labelList& surfacesToTest,
  58     const scalar initDistSqr,
  59     const point& pt
  60 )
  61 {
  62     scalar sum = 0;
  63
  64     forAll(surfacesToTest, testI)
  65     {
  66         label surfI = surfacesToTest[testI];
  67
  68         pointIndexHit hit
  69         (
  70             tempFindNearest(allSurfaces[surfI], pt, initDistSqr)
  71         );
  72
  73         // Note: make it fall over if not hit.
  74         sum += magSqr(hit.hitPoint()-pt);
  75     }
  76     return sum;
  77 }
  78
  79
  80 // Reflects the point furthest away around the triangle centre by a factor fac.
  81 // (triangle centre is the average of all points but the ihi. pSum is running
  82 //  sum of all points)
  83 Foam::scalar Foam::searchableSurfacesQueries::tryMorphTet
  84 (
  85     const PtrList<searchableSurface>& allSurfaces,
  86     const labelList& surfacesToTest,
  87     const scalar initDistSqr,
  88     List<vector>& p,
  89     List<scalar>& y,
  90     vector& pSum,
  91     const label ihi,
  92     const scalar fac
  93 )
  94 {
  95     scalar fac1 = (1.0-fac)/vector::nComponents;
  96     scalar fac2 = fac1-fac;
  97
  98     vector ptry = pSum*fac1-p[ihi]*fac2;
  99
 100     scalar ytry = sumDistSqr(allSurfaces, surfacesToTest, initDistSqr, ptry);
 101
 102     if (ytry < y[ihi])
 103     {
 104         y[ihi] = ytry;
 105         pSum += ptry - p[ihi];
 106         p[ihi] = ptry;
 107     }
 108     return ytry;
 109 }
 110
 111
 112 bool Foam::searchableSurfacesQueries::morphTet
 113 (
 114     const PtrList<searchableSurface>& allSurfaces,
 115     const labelList& surfacesToTest,
 116     const scalar initDistSqr,
 117     const scalar convergenceDistSqr,
 118     const label maxIter,
 119     List<vector>& p,
 120     List<scalar>& y
 121 )
 122 {
 123     vector pSum = sum(p);
 124
 125     autoPtr<OFstream> str;
 126     label vertI = 0;
 127     if (debug)
 128     {
 129         wordList names(surfacesToTest.size());
 130         forAll(surfacesToTest, i)
 131         {
 132             names[i] = allSurfaces[surfacesToTest[i]].name();
 133         }
 134         Pout<< "searchableSurfacesQueries::morphTet : intersection of "
 135             << names << " starting from points:" << p << endl;
 136         str.reset(new OFstream("track.obj"));
 137         meshTools::writeOBJ(str(), p[0]);
 138         vertI++;
 139     }
 140
 141     for (label iter = 0; iter < maxIter; iter++)
 142     {
 143         // Get the indices of lowest, highest and second-highest values.
 144         label ilo, ihi, inhi;
 145         {
 146             labelList sortedIndices;
 147             sortedOrder(y, sortedIndices);
 148             ilo  = sortedIndices[0];
 149             ihi  = sortedIndices[sortedIndices.size()-1];
 150             inhi = sortedIndices[sortedIndices.size()-2];
 151         }
 152
 153         if (debug)
 154         {
 155             Pout<< "Iteration:" << iter
 156                 << " lowest:" << y[ilo] << " highest:" << y[ihi]
 157                 << " points:" << p << endl;
 158
 159             meshTools::writeOBJ(str(), p[ilo]);
 160             vertI++;
 161             str()<< "l " << vertI-1 << ' ' << vertI << nl;
 162         }
 163
 164         if (y[ihi] < convergenceDistSqr)
 165         {
 166             // Get point on 0th surface.
 167             Swap(p[0], p[ilo]);
 168             Swap(y[0], y[ilo]);
 169             return true;
 170         }
 171
 172         // Reflection: point furthest away gets reflected.
 173         scalar ytry = tryMorphTet
 174         (
 175             allSurfaces,
 176             surfacesToTest,
 177             10*y[ihi],             // search box.
 178             p,
 179             y,
 180             pSum,
 181             ihi,
 182             -1.0
 183         );
 184
 185         if (ytry <= y[ilo])
 186         {
 187             // If in right direction (y lower) expand by two.
 188             ytry = tryMorphTet
 189             (
 190                 allSurfaces,
 191                 surfacesToTest,
 192                 10*y[ihi],
 193                 p,
 194                 y,
 195                 pSum,
 196                 ihi,
 197                 2.0
 198             );
 199         }
 200         else if (ytry >= y[inhi])
 201         {
 202             // If inside tet try contraction.
 203
 204             scalar ysave = y[ihi];
 205
 206             ytry = tryMorphTet
 207             (
 208                 allSurfaces,
 209                 surfacesToTest,
 210                 10*y[ihi],
 211                 p,
 212                 y,
 213                 pSum,
 214                 ihi,
 215                 0.5
 216             );
 217
 218             if (ytry >= ysave)
 219             {
 220                 // Contract around lowest point.
 221                 forAll(p, i)
 222                 {
 223                     if (i != ilo)
 224                     {
 225                         p[i] = 0.5*(p[i] + p[ilo]);
 226                         y[i] = sumDistSqr
 227                         (
 228                             allSurfaces,
 229                             surfacesToTest,
 230                             y[ihi],
 231                             p[i]
 232                         );
 233                     }
 234                 }
 235                 pSum = sum(p);
 236             }
 237         }
 238     }
 239
 240     if (debug)
 241     {
 242         meshTools::writeOBJ(str(), p[0]);
 243         vertI++;
 244         str()<< "l " << vertI-1 << ' ' << vertI << nl;
 245     }
 246
 247     // Failure to converge. Return best guess so far.
 248     label ilo = findMin(y);
 249     Swap(p[0], p[ilo]);
 250     Swap(y[0], y[ilo]);
 251     return false;
 252 }
 253
 254
 255 //// Get all intersections (in order) for single surface.
 256 //void Foam::searchableSurfacesQueries::findAllIntersections
 257 //(
 258 //    const searchableSurface& s,
 259 //    const pointField& start,
 260 //    const pointField& end,
 261 //    const vectorField& smallVec,
 262 //    List<List<pointIndexHit> >& surfaceHitInfo
 263 //)
 264 //{
 265 //    surfaceHitInfo.setSize(start.size());
 266 //
 267 //    // Current start point of vector
 268 //    pointField p0(start);
 269 //
 270 //    List<pointIndexHit> intersectInfo(start.size());
 271 //
 272 //    // For test whether finished doing vector.
 273 //    const vectorField dirVec(end-start);
 274 //    const scalarField magSqrDirVec(magSqr(dirVec));
 275 //
 276 //    while (true)
 277 //    {
 278 //        // Find first intersection. Synced.
 279 //        s.findLine(p0, end, intersectInfo);
 280 //
 281 //        label nHits = 0;
 282 //
 283 //        forAll(intersectInfo, i)
 284 //        {
 285 //            if (intersectInfo[i].hit())
 286 //            {
 287 //                nHits++;
 288 //
 289 //                label sz = surfaceHitInfo[i].size();
 290 //                surfaceHitInfo[i].setSize(sz+1);
 291 //                surfaceHitInfo[i][sz] = intersectInfo[i];
 292 //
 293 //                p0[i] = intersectInfo[i].hitPoint() + smallVec[i];
 294 //
 295 //                // If beyond endpoint set to endpoint so as not to pick up
 296 //                // any intersections. Could instead just filter out hits.
 297 //                if (((p0[i]-start[i])&dirVec[i]) > magSqrDirVec[i])
 298 //                {
 299 //                    p0[i] = end[i];
 300 //                }
 301 //            }
 302 //            else
 303 //            {
 304 //                // Set to endpoint to stop intersection test. See above.
 305 //                p0[i] = end[i];
 306 //            }
 307 //        }
 308 //
 309 //        // returnReduce(nHits) ?
 310 //        if (nHits == 0)
 311 //        {
 312 //            break;
 313 //        }
 314 //    }
 315 //}
 316
 317
 318 // Given current set of hits (allSurfaces, allInfo) merge in those coming from
 319 // surface surfI.
 320 void Foam::searchableSurfacesQueries::mergeHits
 321 (
 322     const point& start,
 323     const scalar mergeDist,
 324
 325     const label testI,                          // index of surface
 326     const List<pointIndexHit>& surfHits,  // hits on surface
 327
 328     labelList& allSurfaces,
 329     List<pointIndexHit>& allInfo,
 330     scalarList& allDistSqr
 331 )
 332 {
 333     // Precalculate distances
 334     scalarList surfDistSqr(surfHits.size());
 335     forAll(surfHits, i)
 336     {
 337         surfDistSqr[i] = magSqr(surfHits[i].hitPoint()-start);
 338     }
 339
 340     forAll(surfDistSqr, i)
 341     {
 342         label index = findLower(allDistSqr, surfDistSqr[i]);
 343
 344         // Check if equal to lower.
 345         if
 346         (
 347             index >= 0
 348          && (mag(allDistSqr[index]-surfDistSqr[i]) < mergeDist)
 349         )
 350         {
 351             // Same. Do not count.
 352             //Pout<< "point:" << surfHits[i].hitPoint()
 353             //    << " considered same as:" << allInfo[index].hitPoint()
 354             //    << " within tol:" << mergeDist
 355             //    << endl;
 356         }
 357         else
 358         {
 359             // Check if equal to higher
 360             label next = index+1;
 361             if
 362             (
 363                 next < allDistSqr.size()
 364              && (mag(allDistSqr[next]-surfDistSqr[i]) < mergeDist)
 365             )
 366             {
 367                 //Pout<< "point:" << surfHits[i].hitPoint()
 368                 //    << " considered same as:" << allInfo[next].hitPoint()
 369                 //    << " within tol:" << mergeDist
 370                 //    << endl;
 371             }
 372             else
 373             {
 374                 // Insert after index
 375                 label sz = allSurfaces.size();
 376                 allSurfaces.setSize(sz+1);
 377                 allInfo.setSize(allSurfaces.size());
 378                 allDistSqr.setSize(allSurfaces.size());
 379                 // Make space.
 380                 for (label j = sz-1; j > index; --j)
 381                 {
 382                     allSurfaces[j+1] = allSurfaces[j];
 383                     allInfo[j+1] = allInfo[j];
 384                     allDistSqr[j+1] = allDistSqr[j];
 385                 }
 386                 // Insert new value
 387                 allSurfaces[index+1] = testI;
 388                 allInfo[index+1] = surfHits[i];
 389                 allDistSqr[index+1] = surfDistSqr[i];
 390             }
 391         }
 392     }
 393 }
 394
 395
 396 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 397
 398 // Find any intersection
 399 void Foam::searchableSurfacesQueries::findAnyIntersection
 400 (
 401     const PtrList<searchableSurface>& allSurfaces,
 402     const labelList& surfacesToTest,
 403     const pointField& start,
 404     const pointField& end,
 405     labelList& hitSurfaces,
 406     List<pointIndexHit>& hitInfo
 407 )
 408 {
 409     hitSurfaces.setSize(start.size());
 410     hitSurfaces = -1;
 411     hitInfo.setSize(start.size());
 412
 413     // Work arrays
 414     labelList hitMap(identity(start.size()));
 415     pointField p0(start);
 416     pointField p1(end);
 417     List<pointIndexHit> intersectInfo(start.size());
 418
 419     forAll(surfacesToTest, testI)
 420     {
 421         // Do synchronised call to all surfaces.
 422         allSurfaces[surfacesToTest[testI]].findLineAny(p0, p1, intersectInfo);
 423
 424         // Copy all hits into arguments, continue with misses
 425         label newI = 0;
 426         forAll(intersectInfo, i)
 427         {
 428             if (intersectInfo[i].hit())
 429             {
 430                 hitInfo[hitMap[i]] = intersectInfo[i];
 431                 hitSurfaces[hitMap[i]] = testI;
 432             }
 433             else
 434             {
 435                 if (i != newI)
 436                 {
 437                     hitMap[newI] = hitMap[i];
 438                     p0[newI] = p0[i];
 439                     p1[newI] = p1[i];
 440                 }
 441                 newI++;
 442             }
 443         }
 444
 445         // All done? Note that this decision should be synchronised
 446         if (newI == 0)
 447         {
 448             break;
 449         }
 450
 451         // Trim and continue
 452         hitMap.setSize(newI);
 453         p0.setSize(newI);
 454         p1.setSize(newI);
 455         intersectInfo.setSize(newI);
 456     }
 457 }
 458
 459
 460 void Foam::searchableSurfacesQueries::findAllIntersections
 461 (
 462     const PtrList<searchableSurface>& allSurfaces,
 463     const labelList& surfacesToTest,
 464     const pointField& start,
 465     const pointField& end,
 466     labelListList& hitSurfaces,
 467     List<List<pointIndexHit> >& hitInfo
 468 )
 469 {
 470     // Note: maybe move the single-surface all intersections test into
 471     // searchable surface? Some of the tolerance issues might be
 472     // lessened.
 473
 474     // 2. Currently calling searchableSurface::findLine with start==end
 475     //    is expected to find no intersection. Problem if it does.
 476
 477     hitSurfaces.setSize(start.size());
 478     hitInfo.setSize(start.size());
 479
 480     if (surfacesToTest.empty())
 481     {
 482         return;
 483     }
 484
 485     // Test first surface
 486     allSurfaces[surfacesToTest[0]].findLineAll(start, end, hitInfo);
 487
 488     // Set hitSurfaces and distance
 489     List<scalarList> hitDistSqr(hitInfo.size());
 490     forAll(hitInfo, pointI)
 491     {
 492         const List<pointIndexHit>& pHits = hitInfo[pointI];
 493
 494         labelList& pSurfaces = hitSurfaces[pointI];
 495         pSurfaces.setSize(pHits.size());
 496         pSurfaces = 0;
 497
 498         scalarList& pDistSqr = hitDistSqr[pointI];
 499         pDistSqr.setSize(pHits.size());
 500         forAll(pHits, i)
 501         {
 502             pDistSqr[i] = magSqr(pHits[i].hitPoint() - start[pointI]);
 503         }
 504     }
 505
 506
 507     if (surfacesToTest.size() > 1)
 508     {
 509         // Small vector to increment start vector by to find next intersection
 510         // along line. Constant factor added to make sure that start==end still
 511         // ends iteration in findAllIntersections. Also SMALL is just slightly
 512         // too small.
 513         const vectorField smallVec
 514         (
 515             1E2*SMALL*(end-start)
 516           + vector(ROOTVSMALL,ROOTVSMALL,ROOTVSMALL)
 517         );
 518
 519         // Tolerance used to check whether points are equal. Note: used to
 520         // compare distance^2. Note that we use the maximum possible tolerance
 521         // (reached at intersections close to the end point)
 522         const scalarField mergeDist(2*mag(smallVec)*mag(end-start));
 523
 524         // Test the other surfaces and merge (according to distance from start).
 525         for (label testI = 1; testI < surfacesToTest.size(); testI++)
 526         {
 527             List<List<pointIndexHit> > surfHits;
 528             allSurfaces[surfacesToTest[testI]].findLineAll
 529             (
 530                 start,
 531                 end,
 532                 surfHits
 533             );
 534
 535             forAll(surfHits, pointI)
 536             {
 537                 mergeHits
 538                 (
 539                     start[pointI],          // Current segment
 540                     mergeDist[pointI],
 541
 542                     testI,                  // Surface and its hits
 543                     surfHits[pointI],
 544
 545                     hitSurfaces[pointI],    // Merge into overall hit info
 546                     hitInfo[pointI],
 547                     hitDistSqr[pointI]
 548                 );
 549             }
 550         }
 551     }
 552 }
 553
 554
 555 //// Find intersections of edge nearest to both endpoints.
 556 //void Foam::searchableSurfacesQueries::findNearestIntersection
 557 //(
 558 //    const PtrList<searchableSurface>& allSurfaces,
 559 //    const labelList& surfacesToTest,
 560 //    const pointField& start,
 561 //    const pointField& end,
 562 //
 563 //    labelList& surface1,
 564 //    List<pointIndexHit>& hit1,
 565 //    labelList& surface2,
 566 //    List<pointIndexHit>& hit2
 567 //)
 568 //{
 569 //    // 1. intersection from start to end
 570 //    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 571 //
 572 //    // Initialize arguments
 573 //    surface1.setSize(start.size());
 574 //    surface1 = -1;
 575 //    hit1.setSize(start.size());
 576 //
 577 //    // Current end of segment to test.
 578 //    pointField nearest(end);
 579 //    // Work array
 580 //    List<pointIndexHit> nearestInfo(start.size());
 581 //
 582 //    forAll(surfacesToTest, testI)
 583 //    {
 584 //        // See if any intersection between start and current nearest
 585 //        allSurfaces[surfacesToTest[testI]].findLine
 586 //        (
 587 //            start,
 588 //            nearest,
 589 //            nearestInfo
 590 //        );
 591 //
 592 //        forAll(nearestInfo, pointI)
 593 //        {
 594 //            if (nearestInfo[pointI].hit())
 595 //            {
 596 //                hit1[pointI] = nearestInfo[pointI];
 597 //                surface1[pointI] = testI;
 598 //                nearest[pointI] = hit1[pointI].hitPoint();
 599 //            }
 600 //        }
 601 //    }
 602 //
 603 //
 604 //    // 2. intersection from end to last intersection
 605 //    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 606 //
 607 //    // Find the nearest intersection from end to start. Note that we
 608 //    // initialize to the first intersection (if any).
 609 //    surface2 = surface1;
 610 //    hit2 = hit1;
 611 //
 612 //    // Set current end of segment to test.
 613 //    forAll(nearest, pointI)
 614 //    {
 615 //        if (hit1[pointI].hit())
 616 //        {
 617 //            nearest[pointI] = hit1[pointI].hitPoint();
 618 //        }
 619 //        else
 620 //        {
 621 //            // Disable testing by setting to end.
 622 //            nearest[pointI] = end[pointI];
 623 //        }
 624 //    }
 625 //
 626 //    forAll(surfacesToTest, testI)
 627 //    {
 628 //        // See if any intersection between end and current nearest
 629 //        allSurfaces[surfacesToTest[i]].findLine(end, nearest, nearestInfo);
 630 //
 631 //        forAll(nearestInfo, pointI)
 632 //        {
 633 //            if (nearestInfo[pointI].hit())
 634 //            {
 635 //                hit2[pointI] = nearestInfo[pointI];
 636 //                surface2[pointI] = testI;
 637 //                nearest[pointI] = hit2[pointI].hitPoint();
 638 //            }
 639 //        }
 640 //    }
 641 //}
 642
 643
 644 // Find nearest. Return -1 or nearest point
 645 void Foam::searchableSurfacesQueries::findNearest
 646 (
 647     const PtrList<searchableSurface>& allSurfaces,
 648     const labelList& surfacesToTest,
 649     const pointField& samples,
 650     const scalarField& nearestDistSqr,
 651     labelList& nearestSurfaces,
 652     List<pointIndexHit>& nearestInfo
 653 )
 654 {
 655     // Initialise
 656     nearestSurfaces.setSize(samples.size());
 657     nearestSurfaces = -1;
 658     nearestInfo.setSize(samples.size());
 659
 660     // Work arrays
 661     scalarField minDistSqr(nearestDistSqr);
 662     List<pointIndexHit> hitInfo(samples.size());
 663
 664     forAll(surfacesToTest, testI)
 665     {
 666         allSurfaces[surfacesToTest[testI]].findNearest
 667         (
 668             samples,
 669             minDistSqr,
 670             hitInfo
 671         );
 672
 673         // Update minDistSqr and arguments
 674         forAll(hitInfo, pointI)
 675         {
 676             if (hitInfo[pointI].hit())
 677             {
 678                 minDistSqr[pointI] = magSqr
 679                 (
 680                     hitInfo[pointI].hitPoint()
 681                   - samples[pointI]
 682                 );
 683                 nearestInfo[pointI] = hitInfo[pointI];
 684                 nearestSurfaces[pointI] = testI;
 685             }
 686         }
 687     }
 688 }
 689
 690
 691 //- Calculate point which is on a set of surfaces.
 692 Foam::pointIndexHit Foam::searchableSurfacesQueries::facesIntersection
 693 (
 694     const PtrList<searchableSurface>& allSurfaces,
 695     const labelList& surfacesToTest,
 696     const scalar initDistSqr,
 697     const scalar convergenceDistSqr,
 698     const point& start
 699 )
 700 {
 701     // Get four starting points. Take these as the projection of the
 702     // starting point onto the surfaces and the mid point
 703     List<point> nearest(surfacesToTest.size()+1);
 704
 705     point sumNearest = vector::zero;
 706
 707     forAll(surfacesToTest, i)
 708     {
 709         pointIndexHit hit
 710         (
 711             tempFindNearest(allSurfaces[surfacesToTest[i]], start, initDistSqr)
 712         );
 713
 714         if (hit.hit())
 715         {
 716             nearest[i] = hit.hitPoint();
 717             sumNearest += nearest[i];
 718         }
 719         else
 720         {
 721             FatalErrorIn
 722             (
 723                 "searchableSurfacesQueries::facesIntersection"
 724                 "(const labelList&, const scalar, const scalar, const point&)"
 725             )   << "Did not find point within distance "
 726                 << initDistSqr << " of starting point " << start
 727                 << " on surface "
 728                 << allSurfaces[surfacesToTest[i]].IOobject::name()
 729                 << abort(FatalError);
 730         }
 731     }
 732
 733     nearest[nearest.size()-1] = sumNearest / surfacesToTest.size();
 734
 735
 736     // Get the sum of distances (initial evaluation)
 737     List<scalar> nearestDist(nearest.size());
 738
 739     forAll(nearestDist, i)
 740     {
 741         nearestDist[i] = sumDistSqr
 742         (
 743             allSurfaces,
 744             surfacesToTest,
 745             initDistSqr,
 746             nearest[i]
 747         );
 748     }
 749
 750
 751     //- Downhill Simplex method
 752
 753     bool converged = morphTet
 754     (
 755         allSurfaces,
 756         surfacesToTest,
 757         initDistSqr,
 758         convergenceDistSqr,
 759         2000,
 760         nearest,
 761         nearestDist
 762     );
 763
 764
 765     pointIndexHit intersection;
 766
 767     if (converged)
 768     {
 769         // Project nearest onto 0th surface.
 770         intersection = tempFindNearest
 771         (
 772             allSurfaces[surfacesToTest[0]],
 773             nearest[0],
 774             nearestDist[0]
 775         );
 776     }
 777
 778     //if (!intersection.hit())
 779     //{
 780     //    // Restart
 781     //    scalar smallDist = Foam::sqr(convergenceDistSqr);
 782     //    nearest[0] = intersection.hitPoint();
 783     //    nearest[1] = nearest[0];
 784     //    nearest[1].x() += smallDist;
 785     //    nearest[2] = nearest[0];
 786     //    nearest[2].y() += smallDist;
 787     //    nearest[3] = nearest[0];
 788     //    nearest[3].z() += smallDist;
 789     //
 790     //    forAll(nearestDist, i)
 791     //    {
 792     //        nearestDist[i] = sumDistSqr
 793     //        (
 794     //            surfacesToTest,
 795     //            initDistSqr,
 796     //            nearest[i]
 797     //        );
 798     //    }
 799     //
 800     //    intersection = morphTet
 801     //    (
 802     //        allSurfaces,
 803     //        surfacesToTest,
 804     //        initDistSqr,
 805     //        convergenceDistSqr,
 806     //        1000,
 807     //        nearest,
 808     //        nearestDist
 809     //    );
 810     //}
 811
 812     return intersection;
 813 }
 814
 815
 816
 817 // ************************************************************************* //