src/meshTools/searchableSurface/searchableSurfacesQueries.C

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