src/meshTools/searchableSurface/distributedTriSurfaceMesh.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
   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 "distributedTriSurfaceMesh.H"
  28 #include "mapDistribute.H"
  29 #include "Random.H"
  30 #include "addToRunTimeSelectionTable.H"
  31 #include "triangleFuncs.H"
  32 #include "matchPoints.H"
  33 #include "globalIndex.H"
  34 #include "Time.H"
  35
  36 #include "IFstream.H"
  37 #include "decompositionMethod.H"
  38 #include "vectorList.H"
  39
  40 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  41
  42 namespace Foam
  43 {
  44
  45 defineTypeNameAndDebug(distributedTriSurfaceMesh, 0);
  46 addToRunTimeSelectionTable(searchableSurface, distributedTriSurfaceMesh, dict);
  47
  48 }
  49
  50
  51 template<>
  52 const char*
  53 Foam::NamedEnum<Foam::distributedTriSurfaceMesh::distributionType, 3>::names[] =
  54 {
  55     "follow",
  56     "independent",
  57     "frozen"
  58 };
  59
  60 const Foam::NamedEnum<Foam::distributedTriSurfaceMesh::distributionType, 3>
  61     Foam::distributedTriSurfaceMesh::distributionTypeNames_;
  62
  63
  64 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  65
  66 // Read my additional data from the dictionary
  67 bool Foam::distributedTriSurfaceMesh::read()
  68 {
  69     // Get bb of all domains.
  70     procBb_.setSize(Pstream::nProcs());
  71
  72     procBb_[Pstream::myProcNo()] = List<treeBoundBox>(dict_.lookup("bounds"));
  73     Pstream::gatherList(procBb_);
  74     Pstream::scatterList(procBb_);
  75
  76     // Distribution type
  77     distType_ = distributionTypeNames_.read(dict_.lookup("distributionType"));
  78
  79     // Merge distance
  80     mergeDist_ = readScalar(dict_.lookup("mergeDistance"));
  81
  82     return true;
  83 }
  84
  85
  86 // Is segment fully local?
  87 bool Foam::distributedTriSurfaceMesh::isLocal
  88 (
  89     const List<treeBoundBox>& myBbs,
  90     const point& start,
  91     const point& end
  92 )
  93 {
  94     forAll(myBbs, bbI)
  95     {
  96         if (myBbs[bbI].contains(start) && myBbs[bbI].contains(end))
  97         {
  98             return true;
  99         }
 100     }
 101     return false;
 102 }
 103
 104
 105 //void Foam::distributedTriSurfaceMesh::splitSegment
 106 //(
 107 //    const label segmentI,
 108 //    const point& start,
 109 //    const point& end,
 110 //    const treeBoundBox& bb,
 111 //
 112 //    DynamicList<segment>& allSegments,
 113 //    DynamicList<label>& allSegmentMap,
 114 //    DynamicList<label> sendMap
 115 //) const
 116 //{
 117 //    // Work points
 118 //    point clipPt0, clipPt1;
 119 //
 120 //    if (bb.contains(start))
 121 //    {
 122 //        // start within, trim end to bb
 123 //        bool clipped = bb.intersects(end, start, clipPt0);
 124 //
 125 //        if (clipped)
 126 //        {
 127 //            // segment from start to clippedStart passes
 128 //            // through proc.
 129 //            sendMap[procI].append(allSegments.size());
 130 //            allSegmentMap.append(segmentI);
 131 //            allSegments.append(segment(start, clipPt0));
 132 //        }
 133 //    }
 134 //    else if (bb.contains(end))
 135 //    {
 136 //        // end within, trim start to bb
 137 //        bool clipped = bb.intersects(start, end, clipPt0);
 138 //
 139 //        if (clipped)
 140 //        {
 141 //            sendMap[procI].append(allSegments.size());
 142 //            allSegmentMap.append(segmentI);
 143 //            allSegments.append(segment(clipPt0, end));
 144 //        }
 145 //    }
 146 //    else
 147 //    {
 148 //        // trim both
 149 //        bool clippedStart = bb.intersects(start, end, clipPt0);
 150 //
 151 //        if (clippedStart)
 152 //        {
 153 //            bool clippedEnd = bb.intersects(end, clipPt0, clipPt1);
 154 //
 155 //            if (clippedEnd)
 156 //            {
 157 //                // middle part of segment passes through proc.
 158 //                sendMap[procI].append(allSegments.size());
 159 //                allSegmentMap.append(segmentI);
 160 //                allSegments.append(segment(clipPt0, clipPt1));
 161 //            }
 162 //        }
 163 //    }
 164 //}
 165
 166
 167 void Foam::distributedTriSurfaceMesh::distributeSegment
 168 (
 169     const label segmentI,
 170     const point& start,
 171     const point& end,
 172
 173     DynamicList<segment>& allSegments,
 174     DynamicList<label>& allSegmentMap,
 175     List<DynamicList<label> >& sendMap
 176 ) const
 177 {
 178     // Work points
 179     point clipPt;
 180
 181
 182     // 1. Fully local already handled outside. Note: retest is cheap.
 183     if (isLocal(procBb_[Pstream::myProcNo()], start, end))
 184     {
 185         return;
 186     }
 187
 188
 189     // 2. If fully inside one other processor, then only need to send
 190     // to that one processor even if it intersects another. Rare occurrence
 191     // but cheap to test.
 192     forAll(procBb_, procI)
 193     {
 194         if (procI != Pstream::myProcNo())
 195         {
 196             const List<treeBoundBox>& bbs = procBb_[procI];
 197
 198             if (isLocal(bbs, start, end))
 199             {
 200                 sendMap[procI].append(allSegments.size());
 201                 allSegmentMap.append(segmentI);
 202                 allSegments.append(segment(start, end));
 203                 return;
 204             }
 205         }
 206     }
 207
 208
 209     // 3. If not contained in single processor send to all intersecting
 210     // processors.
 211     forAll(procBb_, procI)
 212     {
 213         const List<treeBoundBox>& bbs = procBb_[procI];
 214
 215         forAll(bbs, bbI)
 216         {
 217             const treeBoundBox& bb = bbs[bbI];
 218
 219             // Scheme a: any processor that intersects the segment gets
 220             // the segment.
 221
 222             if (bb.intersects(start, end, clipPt))
 223             {
 224                 sendMap[procI].append(allSegments.size());
 225                 allSegmentMap.append(segmentI);
 226                 allSegments.append(segment(start, end));
 227             }
 228
 229             // Alternative: any processor only gets clipped bit of
 230             // segment. This gives small problems with additional
 231             // truncation errors.
 232             //splitSegment
 233             //(
 234             //    segmentI,
 235             //    start,
 236             //    end,
 237             //    bb,
 238             //
 239             //    allSegments,
 240             //    allSegmentMap,
 241             //   sendMap[procI]
 242             //);
 243         }
 244     }
 245 }
 246
 247
 248 Foam::autoPtr<Foam::mapDistribute>
 249 Foam::distributedTriSurfaceMesh::distributeSegments
 250 (
 251     const pointField& start,
 252     const pointField& end,
 253
 254     List<segment>& allSegments,
 255     labelList& allSegmentMap
 256 ) const
 257 {
 258     // Determine send map
 259     // ~~~~~~~~~~~~~~~~~~
 260
 261     labelListList sendMap(Pstream::nProcs());
 262
 263     {
 264         // Since intersection test is quite expensive compared to memory
 265         // allocation we use DynamicList to immediately store the segment
 266         // in the correct bin.
 267
 268         // Segments to test
 269         DynamicList<segment> dynAllSegments(start.size());
 270         // Original index of segment
 271         DynamicList<label> dynAllSegmentMap(start.size());
 272         // Per processor indices into allSegments to send
 273         List<DynamicList<label> > dynSendMap(Pstream::nProcs());
 274
 275         forAll(start, segmentI)
 276         {
 277             distributeSegment
 278             (
 279                 segmentI,
 280                 start[segmentI],
 281                 end[segmentI],
 282
 283                 dynAllSegments,
 284                 dynAllSegmentMap,
 285                 dynSendMap
 286             );
 287         }
 288
 289         // Convert dynamicList to labelList
 290         sendMap.setSize(Pstream::nProcs());
 291         forAll(sendMap, procI)
 292         {
 293             dynSendMap[procI].shrink();
 294             sendMap[procI].transfer(dynSendMap[procI]);
 295         }
 296
 297         allSegments.transfer(dynAllSegments.shrink());
 298         allSegmentMap.transfer(dynAllSegmentMap.shrink());
 299     }
 300
 301
 302     // Send over how many I need to receive.
 303     labelListList sendSizes(Pstream::nProcs());
 304     sendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
 305     forAll(sendMap, procI)
 306     {
 307         sendSizes[Pstream::myProcNo()][procI] = sendMap[procI].size();
 308     }
 309     Pstream::gatherList(sendSizes);
 310     Pstream::scatterList(sendSizes);
 311
 312
 313     // Determine order of receiving
 314     labelListList constructMap(Pstream::nProcs());
 315
 316     // My local segments first
 317     constructMap[Pstream::myProcNo()] = identity
 318     (
 319         sendMap[Pstream::myProcNo()].size()
 320     );
 321
 322     label segmentI = constructMap[Pstream::myProcNo()].size();
 323     forAll(constructMap, procI)
 324     {
 325         if (procI != Pstream::myProcNo())
 326         {
 327             // What I need to receive is what other processor is sending to me.
 328             label nRecv = sendSizes[procI][Pstream::myProcNo()];
 329             constructMap[procI].setSize(nRecv);
 330
 331             for (label i = 0; i < nRecv; i++)
 332             {
 333                 constructMap[procI][i] = segmentI++;
 334             }
 335         }
 336     }
 337
 338     return autoPtr<mapDistribute>
 339     (
 340         new mapDistribute
 341         (
 342             segmentI,       // size after construction
 343             sendMap,
 344             constructMap,
 345             true            // reuse storage
 346         )
 347     );
 348 }
 349
 350
 351 void Foam::distributedTriSurfaceMesh::findLine
 352 (
 353     const bool nearestIntersection,
 354     const pointField& start,
 355     const pointField& end,
 356     List<pointIndexHit>& info
 357 ) const
 358 {
 359     const indexedOctree<treeDataTriSurface>& octree = tree();
 360
 361     // Important:force synchronised construction of indexing
 362     const globalIndex& triIndexer = globalTris();
 363
 364     // Initialise
 365     info.setSize(start.size());
 366     forAll(info, i)
 367     {
 368         info[i].setMiss();
 369     }
 370
 371
 372     // Do any local queries
 373     // ~~~~~~~~~~~~~~~~~~~~
 374
 375     label nLocal = 0;
 376
 377     forAll(start, i)
 378     {
 379         if (isLocal(procBb_[Pstream::myProcNo()], start[i], end[i]))
 380         {
 381             if (nearestIntersection)
 382             {
 383                 info[i] = octree.findLine(start[i], end[i]);
 384             }
 385             else
 386             {
 387                 info[i] = octree.findLineAny(start[i], end[i]);
 388             }
 389
 390             if (info[i].hit())
 391             {
 392                 info[i].setIndex(triIndexer.toGlobal(info[i].index()));
 393             }
 394             nLocal++;
 395         }
 396     }
 397
 398
 399     if
 400     (
 401         Pstream::parRun()
 402      && (
 403             returnReduce(nLocal, sumOp<label>())
 404           < returnReduce(start.size(), sumOp<label>())
 405         )
 406     )
 407     {
 408         // Not all can be resolved locally. Build segments and map, send over
 409         // segments, do intersections, send back and merge.
 410
 411
 412         // Construct queries (segments)
 413         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 414
 415         // Segments to test
 416         List<segment> allSegments(start.size());
 417         // Original index of segment
 418         labelList allSegmentMap(start.size());
 419
 420         const autoPtr<mapDistribute> mapPtr
 421         (
 422             distributeSegments
 423             (
 424                 start,
 425                 end,
 426                 allSegments,
 427                 allSegmentMap
 428             )
 429         );
 430         const mapDistribute& map = mapPtr();
 431
 432         label nOldAllSegments = allSegments.size();
 433
 434
 435         // Exchange the segments
 436         // ~~~~~~~~~~~~~~~~~~~~~
 437
 438         map.distribute
 439         (
 440             Pstream::nonBlocking,   //Pstream::scheduled,
 441             List<labelPair>(0),     //map.schedule(),
 442             map.constructSize(),
 443             map.subMap(),           // what to send
 444             map.constructMap(),     // what to receive
 445             allSegments
 446         );
 447
 448
 449         // Do tests I need to do
 450         // ~~~~~~~~~~~~~~~~~~~~~
 451
 452         // Intersections
 453         List<pointIndexHit> intersections(allSegments.size());
 454
 455         forAll(allSegments, i)
 456         {
 457             if (nearestIntersection)
 458             {
 459                 intersections[i] = octree.findLine
 460                 (
 461                     allSegments[i].first(),
 462                     allSegments[i].second()
 463                 );
 464             }
 465             else
 466             {
 467                 intersections[i] = octree.findLineAny
 468                 (
 469                     allSegments[i].first(),
 470                     allSegments[i].second()
 471                 );
 472             }
 473
 474             // Convert triangle index to global numbering
 475             if (intersections[i].hit())
 476             {
 477                 intersections[i].setIndex
 478                 (
 479                     triIndexer.toGlobal(intersections[i].index())
 480                 );
 481             }
 482         }
 483
 484
 485         // Exchange the intersections (opposite to segments)
 486         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 487
 488         map.distribute
 489         (
 490             //Pstream::scheduled,
 491             //map.schedule            // Note reverse schedule
 492             //(
 493             //    map.constructMap(),
 494             //    map.subMap()
 495             //),
 496             Pstream::nonBlocking,
 497             List<labelPair>(0),
 498             nOldAllSegments,
 499             map.constructMap(),     // what to send
 500             map.subMap(),           // what to receive
 501             intersections
 502         );
 503
 504
 505         // Extract the hits
 506         // ~~~~~~~~~~~~~~~~
 507
 508         forAll(intersections, i)
 509         {
 510             const pointIndexHit& allInfo = intersections[i];
 511             label segmentI = allSegmentMap[i];
 512             pointIndexHit& hitInfo = info[segmentI];
 513
 514             if (allInfo.hit())
 515             {
 516                 if (!hitInfo.hit())
 517                 {
 518                     // No intersection yet so take this one
 519                     hitInfo = allInfo;
 520                 }
 521                 else if (nearestIntersection)
 522                 {
 523                     // Nearest intersection
 524                     if
 525                     (
 526                         magSqr(allInfo.hitPoint()-start[segmentI])
 527                       < magSqr(hitInfo.hitPoint()-start[segmentI])
 528                     )
 529                     {
 530                         hitInfo = allInfo;
 531                     }
 532                 }
 533             }
 534         }
 535     }
 536 }
 537
 538
 539 // Exchanges indices to the processor they come from.
 540 // - calculates exchange map
 541 // - uses map to calculate local triangle index
 542 Foam::autoPtr<Foam::mapDistribute>
 543 Foam::distributedTriSurfaceMesh::calcLocalQueries
 544 (
 545     const List<pointIndexHit>& info,
 546     labelList& triangleIndex
 547 ) const
 548 {
 549     triangleIndex.setSize(info.size());
 550
 551     const globalIndex& triIndexer = globalTris();
 552
 553
 554     // Determine send map
 555     // ~~~~~~~~~~~~~~~~~~
 556
 557     // Since determining which processor the query should go to is
 558     // cheap we do a multi-pass algorithm to save some memory temporarily.
 559
 560     // 1. Count
 561     labelList nSend(Pstream::nProcs(), 0);
 562
 563     forAll(info, i)
 564     {
 565         if (info[i].hit())
 566         {
 567             label procI = triIndexer.whichProcID(info[i].index());
 568             nSend[procI]++;
 569         }
 570     }
 571
 572     // 2. Size sendMap
 573     labelListList sendMap(Pstream::nProcs());
 574     forAll(nSend, procI)
 575     {
 576         sendMap[procI].setSize(nSend[procI]);
 577         nSend[procI] = 0;
 578     }
 579
 580     // 3. Fill sendMap
 581     forAll(info, i)
 582     {
 583         if (info[i].hit())
 584         {
 585             label procI = triIndexer.whichProcID(info[i].index());
 586             triangleIndex[i] = triIndexer.toLocal(procI, info[i].index());
 587             sendMap[procI][nSend[procI]++] = i;
 588         }
 589         else
 590         {
 591             triangleIndex[i] = -1;
 592         }
 593     }
 594
 595
 596     // Send over how many I need to receive
 597     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 598
 599     labelListList sendSizes(Pstream::nProcs());
 600     sendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
 601     forAll(sendMap, procI)
 602     {
 603         sendSizes[Pstream::myProcNo()][procI] = sendMap[procI].size();
 604     }
 605     Pstream::gatherList(sendSizes);
 606     Pstream::scatterList(sendSizes);
 607
 608
 609     // Determine receive map
 610     // ~~~~~~~~~~~~~~~~~~~~~
 611
 612     labelListList constructMap(Pstream::nProcs());
 613
 614     // My local segments first
 615     constructMap[Pstream::myProcNo()] = identity
 616     (
 617         sendMap[Pstream::myProcNo()].size()
 618     );
 619
 620     label segmentI = constructMap[Pstream::myProcNo()].size();
 621     forAll(constructMap, procI)
 622     {
 623         if (procI != Pstream::myProcNo())
 624         {
 625             // What I need to receive is what other processor is sending to me.
 626             label nRecv = sendSizes[procI][Pstream::myProcNo()];
 627             constructMap[procI].setSize(nRecv);
 628
 629             for (label i = 0; i < nRecv; i++)
 630             {
 631                 constructMap[procI][i] = segmentI++;
 632             }
 633         }
 634     }
 635
 636
 637     // Pack into distribution map
 638     // ~~~~~~~~~~~~~~~~~~~~~~~~~~
 639
 640     autoPtr<mapDistribute> mapPtr
 641     (
 642         new mapDistribute
 643         (
 644             segmentI,       // size after construction
 645             sendMap,
 646             constructMap,
 647             true            // reuse storage
 648         )
 649     );
 650     const mapDistribute& map = mapPtr();
 651
 652
 653     // Send over queries
 654     // ~~~~~~~~~~~~~~~~~
 655
 656     map.distribute
 657     (
 658         //Pstream::scheduled,
 659         //map.schedule(),
 660         Pstream::nonBlocking,
 661         List<labelPair>(0),
 662         map.constructSize(),
 663         map.subMap(),           // what to send
 664         map.constructMap(),     // what to receive
 665         triangleIndex
 666     );
 667
 668
 669     return mapPtr;
 670 }
 671
 672
 673 Foam::label Foam::distributedTriSurfaceMesh::calcOverlappingProcs
 674 (
 675     const point& centre,
 676     const scalar radiusSqr,
 677     boolList& overlaps
 678 ) const
 679 {
 680     overlaps = false;
 681     label nOverlaps = 0;
 682
 683     forAll(procBb_, procI)
 684     {
 685         const List<treeBoundBox>& bbs = procBb_[procI];
 686
 687         forAll(bbs, bbI)
 688         {
 689             if (bbs[bbI].overlaps(centre, radiusSqr))
 690             {
 691                 overlaps[procI] = true;
 692                 nOverlaps++;
 693                 break;
 694             }
 695         }
 696     }
 697     return nOverlaps;
 698 }
 699
 700
 701 // Generate queries for parallel distance calculation
 702 // - calculates exchange map
 703 // - uses map to exchange points and radius
 704 Foam::autoPtr<Foam::mapDistribute>
 705 Foam::distributedTriSurfaceMesh::calcLocalQueries
 706 (
 707     const pointField& centres,
 708     const scalarField& radiusSqr,
 709
 710     pointField& allCentres,
 711     scalarField& allRadiusSqr,
 712     labelList& allSegmentMap
 713 ) const
 714 {
 715     // Determine queries
 716     // ~~~~~~~~~~~~~~~~~
 717
 718     labelListList sendMap(Pstream::nProcs());
 719
 720     {
 721         // Queries
 722         DynamicList<point> dynAllCentres(centres.size());
 723         DynamicList<scalar> dynAllRadiusSqr(centres.size());
 724         // Original index of segment
 725         DynamicList<label> dynAllSegmentMap(centres.size());
 726         // Per processor indices into allSegments to send
 727         List<DynamicList<label> > dynSendMap(Pstream::nProcs());
 728
 729         // Work array - whether processor bb overlaps the bounding sphere.
 730         boolList procBbOverlaps(Pstream::nProcs());
 731
 732         forAll(centres, centreI)
 733         {
 734             // Find the processor this sample+radius overlaps.
 735             calcOverlappingProcs
 736             (
 737                 centres[centreI],
 738                 radiusSqr[centreI],
 739                 procBbOverlaps
 740             );
 741
 742             forAll(procBbOverlaps, procI)
 743             {
 744                 if (procI != Pstream::myProcNo() && procBbOverlaps[procI])
 745                 {
 746                     dynSendMap[procI].append(dynAllCentres.size());
 747                     dynAllSegmentMap.append(centreI);
 748                     dynAllCentres.append(centres[centreI]);
 749                     dynAllRadiusSqr.append(radiusSqr[centreI]);
 750                 }
 751             }
 752         }
 753
 754         // Convert dynamicList to labelList
 755         sendMap.setSize(Pstream::nProcs());
 756         forAll(sendMap, procI)
 757         {
 758             dynSendMap[procI].shrink();
 759             sendMap[procI].transfer(dynSendMap[procI]);
 760         }
 761
 762         allCentres.transfer(dynAllCentres.shrink());
 763         allRadiusSqr.transfer(dynAllRadiusSqr.shrink());
 764         allSegmentMap.transfer(dynAllSegmentMap.shrink());
 765     }
 766
 767
 768     // Send over how many I need to receive.
 769     labelListList sendSizes(Pstream::nProcs());
 770     sendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
 771     forAll(sendMap, procI)
 772     {
 773         sendSizes[Pstream::myProcNo()][procI] = sendMap[procI].size();
 774     }
 775     Pstream::gatherList(sendSizes);
 776     Pstream::scatterList(sendSizes);
 777
 778
 779     // Determine order of receiving
 780     labelListList constructMap(Pstream::nProcs());
 781
 782     // My local segments first
 783     constructMap[Pstream::myProcNo()] = identity
 784     (
 785         sendMap[Pstream::myProcNo()].size()
 786     );
 787
 788     label segmentI = constructMap[Pstream::myProcNo()].size();
 789     forAll(constructMap, procI)
 790     {
 791         if (procI != Pstream::myProcNo())
 792         {
 793             // What I need to receive is what other processor is sending to me.
 794             label nRecv = sendSizes[procI][Pstream::myProcNo()];
 795             constructMap[procI].setSize(nRecv);
 796
 797             for (label i = 0; i < nRecv; i++)
 798             {
 799                 constructMap[procI][i] = segmentI++;
 800             }
 801         }
 802     }
 803
 804     autoPtr<mapDistribute> mapPtr
 805     (
 806         new mapDistribute
 807         (
 808             segmentI,       // size after construction
 809             sendMap,
 810             constructMap,
 811             true            // reuse storage
 812         )
 813     );
 814     return mapPtr;
 815 }
 816
 817
 818 // Find bounding boxes that guarantee a more or less uniform distribution
 819 // of triangles. Decomposition in here is only used to get the bounding
 820 // boxes, actual decomposition is done later on.
 821 // Returns a per processor a list of bounding boxes that most accurately
 822 // describe the shape. For now just a single bounding box per processor but
 823 // optimisation might be to determine a better fitting shape.
 824 Foam::List<Foam::List<Foam::treeBoundBox> >
 825 Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
 826 (
 827     const triSurface& s
 828 )
 829 {
 830     if (!decomposer_.valid())
 831     {
 832         // Use current decomposer.
 833         // Note: or always use hierarchical?
 834         IOdictionary decomposeDict
 835         (
 836             IOobject
 837             (
 838                 "decomposeParDict",
 839                 searchableSurface::time().system(),
 840                 searchableSurface::time(),
 841                 IOobject::MUST_READ,
 842                 IOobject::NO_WRITE,
 843                 false
 844             )
 845         );
 846         decomposer_ = decompositionMethod::New(decomposeDict);
 847
 848         if (!decomposer_().parallelAware())
 849         {
 850             FatalErrorIn
 851             (
 852                 "distributedTriSurfaceMesh::independentlyDistributedBbs"
 853                 "(const triSurface&)"
 854             )   << "The decomposition method " << decomposer_().typeName
 855                 << " does not decompose in parallel."
 856                 << " Please choose one that does." << exit(FatalError);
 857         }
 858     }
 859
 860     // Do decomposition according to triangle centre
 861     pointField triCentres(s.size());
 862     forAll (s, triI)
 863     {
 864         triCentres[triI] = s[triI].centre(s.points());
 865     }
 866
 867     // Do the actual decomposition
 868     labelList distribution(decomposer_->decompose(triCentres));
 869
 870     // Find bounding box for all triangles on new distribution.
 871
 872     // Initialise to inverted box (VGREAT, -VGREAT)
 873     List<List<treeBoundBox> > bbs(Pstream::nProcs());
 874     forAll(bbs, procI)
 875     {
 876         bbs[procI].setSize(1);
 877         //bbs[procI][0] = boundBox::invertedBox;
 878         bbs[procI][0].min() = point( VGREAT,  VGREAT,  VGREAT);
 879         bbs[procI][0].max() = point(-VGREAT, -VGREAT, -VGREAT);
 880     }
 881
 882     forAll (s, triI)
 883     {
 884         point& bbMin = bbs[distribution[triI]][0].min();
 885         point& bbMax = bbs[distribution[triI]][0].max();
 886
 887         const labelledTri& f = s[triI];
 888         const point& p0 = s.points()[f[0]];
 889         const point& p1 = s.points()[f[1]];
 890         const point& p2 = s.points()[f[2]];
 891
 892         bbMin = min(bbMin, p0);
 893         bbMin = min(bbMin, p1);
 894         bbMin = min(bbMin, p2);
 895
 896         bbMax = max(bbMax, p0);
 897         bbMax = max(bbMax, p1);
 898         bbMax = max(bbMax, p2);
 899     }
 900
 901     // Now combine for all processors and convert to correct format.
 902     forAll(bbs, procI)
 903     {
 904         forAll(bbs[procI], i)
 905         {
 906             reduce(bbs[procI][i].min(), minOp<point>());
 907             reduce(bbs[procI][i].max(), maxOp<point>());
 908         }
 909     }
 910     return bbs;
 911 }
 912
 913
 914 void Foam::distributedTriSurfaceMesh::calcBounds
 915 (
 916     boundBox& bb,
 917     label& nPoints
 918 ) const
 919 {
 920     // Unfortunately nPoints constructs meshPoints() so do compact version
 921     // ourselves
 922
 923     PackedList<1> pointIsUsed(points().size());
 924     pointIsUsed = 0U;
 925
 926     nPoints = 0;
 927     bb.min() = point(VGREAT, VGREAT, VGREAT);
 928     bb.max() = point(-VGREAT, -VGREAT, -VGREAT);
 929
 930     const triSurface& s = static_cast<const triSurface&>(*this);
 931
 932     forAll(s, triI)
 933     {
 934         const labelledTri& f = s[triI];
 935
 936         forAll(f, fp)
 937         {
 938             label pointI = f[fp];
 939             if (pointIsUsed.set(pointI, 1))
 940             {
 941                 bb.min() = ::Foam::min(bb.min(), points()[pointI]);
 942                 bb.max() = ::Foam::max(bb.max(), points()[pointI]);
 943                 nPoints++;
 944             }
 945         }
 946     }
 947 }
 948
 949
 950 // Does any part of triangle overlap bb.
 951 bool Foam::distributedTriSurfaceMesh::overlaps
 952 (
 953     const List<treeBoundBox>& bbs,
 954     const point& p0,
 955     const point& p1,
 956     const point& p2
 957 )
 958 {
 959     forAll(bbs, bbI)
 960     {
 961         const treeBoundBox& bb = bbs[bbI];
 962
 963         boundBox triBb(p0, p0);
 964         triBb.min() = min(triBb.min(), p1);
 965         triBb.min() = min(triBb.min(), p2);
 966
 967         triBb.max() = max(triBb.max(), p1);
 968         triBb.max() = max(triBb.max(), p2);
 969
 970         //- Exact test of triangle intersecting bb
 971
 972         // Quick rejection. If whole bounding box of tri is outside cubeBb then
 973         // there will be no intersection.
 974         if (bb.overlaps(triBb))
 975         {
 976             // Check if one or more triangle point inside
 977             if (bb.contains(p0) || bb.contains(p1) || bb.contains(p2))
 978             {
 979                 // One or more points inside
 980                 return true;
 981             }
 982
 983             // Now we have the difficult case: all points are outside but
 984             // connecting edges might go through cube. Use fast intersection
 985             // of bounding box.
 986             bool intersect = triangleFuncs::intersectBb(p0, p1, p2, bb);
 987
 988             if (intersect)
 989             {
 990                 return true;
 991             }
 992         }
 993     }
 994     return false;
 995 }
 996
 997
 998 void Foam::distributedTriSurfaceMesh::subsetMeshMap
 999 (
1000     const triSurface& s,
1001     const boolList& include,
1002     const label nIncluded,
1003     labelList& newToOldPoints,
1004     labelList& oldToNewPoints,
1005     labelList& newToOldFaces
1006 )
1007 {
1008     newToOldFaces.setSize(nIncluded);
1009     newToOldPoints.setSize(s.points().size());
1010     oldToNewPoints.setSize(s.points().size());
1011     oldToNewPoints = -1;
1012     {
1013         label faceI = 0;
1014         label pointI = 0;
1015
1016         forAll(include, oldFacei)
1017         {
1018             if (include[oldFacei])
1019             {
1020                 // Store new faces compact
1021                 newToOldFaces[faceI++] = oldFacei;
1022
1023                 // Renumber labels for triangle
1024                 const labelledTri& tri = s[oldFacei];
1025
1026                 forAll(tri, fp)
1027                 {
1028                     label oldPointI = tri[fp];
1029
1030                     if (oldToNewPoints[oldPointI] == -1)
1031                     {
1032                         oldToNewPoints[oldPointI] = pointI;
1033                         newToOldPoints[pointI++] = oldPointI;
1034                     }
1035                 }
1036             }
1037         }
1038         newToOldPoints.setSize(pointI);
1039     }
1040 }
1041
1042
1043 Foam::triSurface Foam::distributedTriSurfaceMesh::subsetMesh
1044 (
1045     const triSurface& s,
1046     const labelList& newToOldPoints,
1047     const labelList& oldToNewPoints,
1048     const labelList& newToOldFaces
1049 )
1050 {
1051     // Extract points
1052     pointField newPoints(newToOldPoints.size());
1053     forAll(newToOldPoints, i)
1054     {
1055         newPoints[i] = s.points()[newToOldPoints[i]];
1056     }
1057     // Extract faces
1058     List<labelledTri> newTriangles(newToOldFaces.size());
1059
1060     forAll(newToOldFaces, i)
1061     {
1062         // Get old vertex labels
1063         const labelledTri& tri = s[newToOldFaces[i]];
1064
1065         newTriangles[i][0] = oldToNewPoints[tri[0]];
1066         newTriangles[i][1] = oldToNewPoints[tri[1]];
1067         newTriangles[i][2] = oldToNewPoints[tri[2]];
1068         newTriangles[i].region() = tri.region();
1069     }
1070
1071     // Reuse storage.
1072     return triSurface(newTriangles, s.patches(), newPoints, true);
1073 }
1074
1075
1076 Foam::triSurface Foam::distributedTriSurfaceMesh::subsetMesh
1077 (
1078     const triSurface& s,
1079     const boolList& include,
1080     labelList& newToOldPoints,
1081     labelList& newToOldFaces
1082 )
1083 {
1084     label n = 0;
1085
1086     forAll(include, i)
1087     {
1088         if (include[i])
1089         {
1090             n++;
1091         }
1092     }
1093
1094     labelList oldToNewPoints;
1095     subsetMeshMap
1096     (
1097         s,
1098         include,
1099         n,
1100         newToOldPoints,
1101         oldToNewPoints,
1102         newToOldFaces
1103     );
1104
1105     return subsetMesh
1106     (
1107         s,
1108         newToOldPoints,
1109         oldToNewPoints,
1110         newToOldFaces
1111     );
1112 }
1113
1114
1115 Foam::triSurface Foam::distributedTriSurfaceMesh::subsetMesh
1116 (
1117     const triSurface& s,
1118     const labelList& newToOldFaces,
1119     labelList& newToOldPoints
1120 )
1121 {
1122     const boolList include
1123     (
1124         createWithValues<boolList>
1125         (
1126             s.size(),
1127             false,
1128             newToOldFaces,
1129             true
1130         )
1131     );
1132
1133     newToOldPoints.setSize(s.points().size());
1134     labelList oldToNewPoints(s.points().size(), -1);
1135     {
1136         label pointI = 0;
1137
1138         forAll(include, oldFacei)
1139         {
1140             if (include[oldFacei])
1141             {
1142                 // Renumber labels for triangle
1143                 const labelledTri& tri = s[oldFacei];
1144
1145                 forAll(tri, fp)
1146                 {
1147                     label oldPointI = tri[fp];
1148
1149                     if (oldToNewPoints[oldPointI] == -1)
1150                     {
1151                         oldToNewPoints[oldPointI] = pointI;
1152                         newToOldPoints[pointI++] = oldPointI;
1153                     }
1154                 }
1155             }
1156         }
1157         newToOldPoints.setSize(pointI);
1158     }
1159
1160     return subsetMesh
1161     (
1162         s,
1163         newToOldPoints,
1164         oldToNewPoints,
1165         newToOldFaces
1166     );
1167 }
1168
1169
1170 Foam::label Foam::distributedTriSurfaceMesh::findTriangle
1171 (
1172     const List<labelledTri>& allFaces,
1173     const labelListList& allPointFaces,
1174     const labelledTri& otherF
1175 )
1176 {
1177     // allFaces connected to otherF[0]
1178     const labelList& pFaces = allPointFaces[otherF[0]];
1179
1180     forAll(pFaces, i)
1181     {
1182         const labelledTri& f = allFaces[pFaces[i]];
1183
1184         if (f.region() == otherF.region())
1185         {
1186             // Find index of otherF[0]
1187             label fp0 = findIndex(f, otherF[0]);
1188             // Check rest of triangle in same order
1189             label fp1 = f.fcIndex(fp0);
1190             label fp2 = f.fcIndex(fp1);
1191
1192             if (f[fp1] == otherF[1] && f[fp2] == otherF[2])
1193             {
1194                 return pFaces[i];
1195             }
1196         }
1197     }
1198     return -1;
1199 }
1200
1201
1202 // Merge into allSurf.
1203 void Foam::distributedTriSurfaceMesh::merge
1204 (
1205     const scalar mergeDist,
1206     const List<labelledTri>& subTris,
1207     const pointField& subPoints,
1208
1209     List<labelledTri>& allTris,
1210     pointField& allPoints,
1211
1212     labelList& faceConstructMap,
1213     labelList& pointConstructMap
1214 )
1215 {
1216     labelList subToAll;
1217     matchPoints
1218     (
1219         subPoints,
1220         allPoints,
1221         scalarField(subPoints.size(), mergeDist),   // match distance
1222         false,                                      // verbose
1223         pointConstructMap
1224     );
1225
1226     label nOldAllPoints = allPoints.size();
1227
1228
1229     // Add all unmatched points
1230     // ~~~~~~~~~~~~~~~~~~~~~~~~
1231
1232     label allPointI = nOldAllPoints;
1233     forAll(pointConstructMap, pointI)
1234     {
1235         if (pointConstructMap[pointI] == -1)
1236         {
1237             pointConstructMap[pointI] = allPointI++;
1238         }
1239     }
1240
1241     if (allPointI > nOldAllPoints)
1242     {
1243         allPoints.setSize(allPointI);
1244
1245         forAll(pointConstructMap, pointI)
1246         {
1247             if (pointConstructMap[pointI] >= nOldAllPoints)
1248             {
1249                 allPoints[pointConstructMap[pointI]] = subPoints[pointI];
1250             }
1251         }
1252     }
1253
1254
1255     // To check whether triangles are same we use pointFaces.
1256     labelListList allPointFaces;
1257     invertManyToMany(nOldAllPoints, allTris, allPointFaces);
1258
1259
1260     // Add all unmatched triangles
1261     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
1262
1263     label allTriI = allTris.size();
1264     allTris.setSize(allTriI + subTris.size());
1265
1266     faceConstructMap.setSize(subTris.size());
1267
1268     forAll(subTris, triI)
1269     {
1270         const labelledTri& subTri = subTris[triI];
1271
1272         // Get triangle in new numbering
1273         labelledTri mappedTri
1274         (
1275             pointConstructMap[subTri[0]],
1276             pointConstructMap[subTri[1]],
1277             pointConstructMap[subTri[2]],
1278             subTri.region()
1279         );
1280
1281
1282         // Check if all points were already in surface
1283         bool fullMatch = true;
1284
1285         forAll(mappedTri, fp)
1286         {
1287             if (mappedTri[fp] >= nOldAllPoints)
1288             {
1289                 fullMatch = false;
1290                 break;
1291             }
1292         }
1293
1294         if (fullMatch)
1295         {
1296             // All three points are mapped to old points. See if same
1297             // triangle.
1298             label i = findTriangle
1299             (
1300                 allTris,
1301                 allPointFaces,
1302                 mappedTri
1303             );
1304
1305             if (i == -1)
1306             {
1307                 // Add
1308                 faceConstructMap[triI] = allTriI;
1309                 allTris[allTriI] = mappedTri;
1310                 allTriI++;
1311             }
1312             else
1313             {
1314                 faceConstructMap[triI] = i;
1315             }
1316         }
1317         else
1318         {
1319             // Add
1320             faceConstructMap[triI] = allTriI;
1321             allTris[allTriI] = mappedTri;
1322             allTriI++;
1323         }
1324     }
1325     allTris.setSize(allTriI);
1326 }
1327
1328
1329 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
1330
1331 Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh
1332 (
1333     const IOobject& io,
1334     const triSurface& s,
1335     const dictionary& dict
1336 )
1337 :
1338     triSurfaceMesh(io, s),
1339     dict_
1340     (
1341         IOobject
1342         (
1343             searchableSurface::name() + "Dict",
1344             searchableSurface::instance(),
1345             searchableSurface::local(),
1346             searchableSurface::db(),
1347             searchableSurface::NO_READ,
1348             searchableSurface::writeOpt(),
1349             searchableSurface::registerObject()
1350         ),
1351         dict
1352     )
1353 {
1354     read();
1355
1356     if (debug)
1357     {
1358         Info<< "Constructed from triSurface:" << endl;
1359         writeStats(Info);
1360
1361         labelList nTris(Pstream::nProcs());
1362         nTris[Pstream::myProcNo()] = triSurface::size();
1363         Pstream::gatherList(nTris);
1364         Pstream::scatterList(nTris);
1365
1366         Info<< endl<< "\tproc\ttris\tbb" << endl;
1367         forAll(nTris, procI)
1368         {
1369             Info<< '\t' << procI << '\t' << nTris[procI]
1370                  << '\t' << procBb_[procI] << endl;
1371         }
1372         Info<< endl;
1373     }
1374 }
1375
1376
1377 Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh(const IOobject& io)
1378 :
1379     //triSurfaceMesh(io),
1380     triSurfaceMesh
1381     (
1382         IOobject
1383         (
1384             io.name(),
1385             io.time().findInstance(io.local(), word::null),
1386             io.local(),
1387             io.db(),
1388             io.readOpt(),
1389             io.writeOpt(),
1390             io.registerObject()
1391         )
1392     ),
1393     dict_
1394     (
1395         IOobject
1396         (
1397             searchableSurface::name() + "Dict",
1398             searchableSurface::instance(),
1399             searchableSurface::local(),
1400             searchableSurface::db(),
1401             searchableSurface::readOpt(),
1402             searchableSurface::writeOpt(),
1403             searchableSurface::registerObject()
1404         )
1405     )
1406 {
1407     read();
1408
1409     if (debug)
1410     {
1411         Info<< "Read distributedTriSurface from " << io.objectPath()
1412             << ':' << endl;
1413         writeStats(Info);
1414
1415         labelList nTris(Pstream::nProcs());
1416         nTris[Pstream::myProcNo()] = triSurface::size();
1417         Pstream::gatherList(nTris);
1418         Pstream::scatterList(nTris);
1419
1420         Info<< endl<< "\tproc\ttris\tbb" << endl;
1421         forAll(nTris, procI)
1422         {
1423             Info<< '\t' << procI << '\t' << nTris[procI]
1424                  << '\t' << procBb_[procI] << endl;
1425         }
1426         Info<< endl;
1427     }
1428 }
1429
1430
1431 Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh
1432 (
1433     const IOobject& io,
1434     const dictionary& dict
1435 )
1436 :
1437     //triSurfaceMesh(io, dict),
1438     triSurfaceMesh
1439     (
1440         IOobject
1441         (
1442             io.name(),
1443             io.time().findInstance(io.local(), word::null),
1444             io.local(),
1445             io.db(),
1446             io.readOpt(),
1447             io.writeOpt(),
1448             io.registerObject()
1449         ),
1450         dict
1451     ),
1452     dict_
1453     (
1454         IOobject
1455         (
1456             searchableSurface::name() + "Dict",
1457             searchableSurface::instance(),
1458             searchableSurface::local(),
1459             searchableSurface::db(),
1460             searchableSurface::readOpt(),
1461             searchableSurface::writeOpt(),
1462             searchableSurface::registerObject()
1463         )
1464     )
1465 {
1466     read();
1467
1468     if (debug)
1469     {
1470         Info<< "Read distributedTriSurface from " << io.objectPath()
1471             << " and dictionary:" << endl;
1472         writeStats(Info);
1473
1474         labelList nTris(Pstream::nProcs());
1475         nTris[Pstream::myProcNo()] = triSurface::size();
1476         Pstream::gatherList(nTris);
1477         Pstream::scatterList(nTris);
1478
1479         Info<< endl<< "\tproc\ttris\tbb" << endl;
1480         forAll(nTris, procI)
1481         {
1482             Info<< '\t' << procI << '\t' << nTris[procI]
1483                  << '\t' << procBb_[procI] << endl;
1484         }
1485         Info<< endl;
1486     }
1487 }
1488
1489
1490 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
1491
1492 Foam::distributedTriSurfaceMesh::~distributedTriSurfaceMesh()
1493 {
1494     clearOut();
1495 }
1496
1497
1498 void Foam::distributedTriSurfaceMesh::clearOut()
1499 {
1500     globalTris_.clear();
1501     triSurfaceMesh::clearOut();
1502 }
1503
1504
1505 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
1506
1507 const Foam::globalIndex& Foam::distributedTriSurfaceMesh::globalTris() const
1508 {
1509     if (!globalTris_.valid())
1510     {
1511         globalTris_.reset(new globalIndex(triSurface::size()));
1512     }
1513     return globalTris_;
1514 }
1515
1516
1517 void Foam::distributedTriSurfaceMesh::findNearest
1518 (
1519     const pointField& samples,
1520     const scalarField& nearestDistSqr,
1521     List<pointIndexHit>& info
1522 ) const
1523 {
1524     const indexedOctree<treeDataTriSurface>& octree = tree();
1525
1526     // Important:force synchronised construction of indexing
1527     const globalIndex& triIndexer = globalTris();
1528
1529
1530     // Initialise
1531     // ~~~~~~~~~~
1532
1533     info.setSize(samples.size());
1534     forAll(info, i)
1535     {
1536         info[i].setMiss();
1537     }
1538
1539
1540
1541     // Do any local queries
1542     // ~~~~~~~~~~~~~~~~~~~~
1543
1544     label nLocal = 0;
1545
1546     {
1547         // Work array - whether processor bb overlaps the bounding sphere.
1548         boolList procBbOverlaps(Pstream::nProcs());
1549
1550         forAll(samples, i)
1551         {
1552             // Find the processor this sample+radius overlaps.
1553             label nProcs = calcOverlappingProcs
1554             (
1555                 samples[i],
1556                 nearestDistSqr[i],
1557                 procBbOverlaps
1558             );
1559
1560             // Overlaps local processor?
1561             if (procBbOverlaps[Pstream::myProcNo()])
1562             {
1563                 info[i] = octree.findNearest(samples[i], nearestDistSqr[i]);
1564                 if (info[i].hit())
1565                 {
1566                     info[i].setIndex(triIndexer.toGlobal(info[i].index()));
1567                 }
1568                 if (nProcs == 1)
1569                 {
1570                     // Fully local
1571                     nLocal++;
1572                 }
1573             }
1574         }
1575     }
1576
1577
1578     if
1579     (
1580         Pstream::parRun()
1581      && (
1582             returnReduce(nLocal, sumOp<label>())
1583           < returnReduce(samples.size(), sumOp<label>())
1584         )
1585     )
1586     {
1587         // Not all can be resolved locally. Build queries and map, send over
1588         // queries, do intersections, send back and merge.
1589
1590         // Calculate queries and exchange map
1591         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1592
1593         pointField allCentres;
1594         scalarField allRadiusSqr;
1595         labelList allSegmentMap;
1596         autoPtr<mapDistribute> mapPtr
1597         (
1598             calcLocalQueries
1599             (
1600                 samples,
1601                 nearestDistSqr,
1602
1603                 allCentres,
1604                 allRadiusSqr,
1605                 allSegmentMap
1606             )
1607         );
1608         const mapDistribute& map = mapPtr();
1609
1610
1611         // swap samples to local processor
1612         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1613
1614         map.distribute
1615         (
1616             //Pstream::scheduled,
1617             //map.schedule(),
1618             Pstream::nonBlocking,
1619             List<labelPair>(0),
1620             map.constructSize(),
1621             map.subMap(),           // what to send
1622             map.constructMap(),     // what to receive
1623             allCentres
1624         );
1625         map.distribute
1626         (
1627             //Pstream::scheduled,
1628             //map.schedule(),
1629             Pstream::nonBlocking,
1630             List<labelPair>(0),
1631             map.constructSize(),
1632             map.subMap(),           // what to send
1633             map.constructMap(),     // what to receive
1634             allRadiusSqr
1635         );
1636
1637
1638         // Do my tests
1639         // ~~~~~~~~~~~
1640
1641         List<pointIndexHit> allInfo(allCentres.size());
1642         forAll(allInfo, i)
1643         {
1644             allInfo[i] = octree.findNearest
1645             (
1646                 allCentres[i],
1647                 allRadiusSqr[i]
1648             );
1649             if (allInfo[i].hit())
1650             {
1651                 allInfo[i].setIndex(triIndexer.toGlobal(allInfo[i].index()));
1652             }
1653         }
1654
1655
1656         // Send back results
1657         // ~~~~~~~~~~~~~~~~~
1658
1659         map.distribute
1660         (
1661             //Pstream::scheduled,
1662             //map.schedule            // note reverse schedule
1663             //(
1664             //    map.constructMap(),
1665             //    map.subMap()
1666             //),
1667             Pstream::nonBlocking,
1668             List<labelPair>(0),
1669             allSegmentMap.size(),
1670             map.constructMap(),     // what to send
1671             map.subMap(),           // what to receive
1672             allInfo
1673         );
1674
1675
1676         // Extract information
1677         // ~~~~~~~~~~~~~~~~~~~
1678
1679         forAll(allInfo, i)
1680         {
1681             if (allInfo[i].hit())
1682             {
1683                 label pointI = allSegmentMap[i];
1684
1685                 if (!info[pointI].hit())
1686                 {
1687                     // No intersection yet so take this one
1688                     info[pointI] = allInfo[i];
1689                 }
1690                 else
1691                 {
1692                     // Nearest intersection
1693                     if
1694                     (
1695                         magSqr(allInfo[i].hitPoint()-samples[pointI])
1696                       < magSqr(info[pointI].hitPoint()-samples[pointI])
1697                     )
1698                     {
1699                         info[pointI] = allInfo[i];
1700                     }
1701                 }
1702             }
1703         }
1704     }
1705 }
1706
1707
1708 void Foam::distributedTriSurfaceMesh::findLine
1709 (
1710     const pointField& start,
1711     const pointField& end,
1712     List<pointIndexHit>& info
1713 ) const
1714 {
1715     findLine
1716     (
1717         true,   // nearestIntersection
1718         start,
1719         end,
1720         info
1721     );
1722 }
1723
1724
1725 void Foam::distributedTriSurfaceMesh::findLineAny
1726 (
1727     const pointField& start,
1728     const pointField& end,
1729     List<pointIndexHit>& info
1730 ) const
1731 {
1732     findLine
1733     (
1734         true,   // nearestIntersection
1735         start,
1736         end,
1737         info
1738     );
1739 }
1740
1741
1742 void Foam::distributedTriSurfaceMesh::findLineAll
1743 (
1744     const pointField& start,
1745     const pointField& end,
1746     List<List<pointIndexHit> >& info
1747 ) const
1748 {
1749     // Reuse fineLine. We could modify all of findLine to do multiple
1750     // intersections but this would mean a lot of data transferred so
1751     // for now we just find nearest intersection and retest from that
1752     // intersection onwards.
1753
1754     // Work array.
1755     List<pointIndexHit> hitInfo(start.size());
1756
1757     findLine
1758     (
1759         true,   // nearestIntersection
1760         start,
1761         end,
1762         hitInfo
1763     );
1764
1765     // Tolerances:
1766     // To find all intersections we add a small vector to the last intersection
1767     // This is chosen such that
1768     // - it is significant (SMALL is smallest representative relative tolerance;
1769     //   we need something bigger since we're doing calculations)
1770     // - if the start-end vector is zero we still progress
1771     const vectorField dirVec(end-start);
1772     const scalarField magSqrDirVec(magSqr(dirVec));
1773     const vectorField smallVec
1774     (
1775         Foam::sqrt(SMALL)*dirVec
1776       + vector(ROOTVSMALL,ROOTVSMALL,ROOTVSMALL)
1777     );
1778
1779     // Copy to input and compact any hits
1780     labelList pointMap(start.size());
1781     pointField e0(start.size());
1782     pointField e1(start.size());
1783     label compactI = 0;
1784
1785     info.setSize(hitInfo.size());
1786     forAll(hitInfo, pointI)
1787     {
1788         if (hitInfo[pointI].hit())
1789         {
1790             info[pointI].setSize(1);
1791             info[pointI][0] = hitInfo[pointI];
1792
1793             point pt = hitInfo[pointI].hitPoint() + smallVec[pointI];
1794
1795             if (((pt-start[pointI])&dirVec[pointI]) <= magSqrDirVec[pointI])
1796             {
1797                 e0[compactI] = pt;
1798                 e1[compactI] = end[pointI];
1799                 pointMap[compactI] = pointI;
1800                 compactI++;
1801             }
1802         }
1803         else
1804         {
1805             info[pointI].clear();
1806         }
1807     }
1808
1809     e0.setSize(compactI);
1810     e1.setSize(compactI);
1811     pointMap.setSize(compactI);
1812
1813     while (returnReduce(e0.size(), sumOp<label>()) > 0)
1814     {
1815         findLine
1816         (
1817             true,   // nearestIntersection
1818             e0,
1819             e1,
1820             hitInfo
1821         );
1822
1823
1824         // Extract
1825         label compactI = 0;
1826         forAll(hitInfo, i)
1827         {
1828             if (hitInfo[i].hit())
1829             {
1830                 label pointI = pointMap[i];
1831
1832                 label sz = info[pointI].size();
1833                 info[pointI].setSize(sz+1);
1834                 info[pointI][sz] = hitInfo[i];
1835
1836                 point pt = hitInfo[i].hitPoint() + smallVec[pointI];
1837
1838                 if (((pt-start[pointI])&dirVec[pointI]) <= magSqrDirVec[pointI])
1839                 {
1840                     e0[compactI] = pt;
1841                     e1[compactI] = end[pointI];
1842                     pointMap[compactI] = pointI;
1843                     compactI++;
1844                 }
1845             }
1846         }
1847
1848         // Trim
1849         e0.setSize(compactI);
1850         e1.setSize(compactI);
1851         pointMap.setSize(compactI);
1852     }
1853 }
1854
1855
1856 void Foam::distributedTriSurfaceMesh::getRegion
1857 (
1858     const List<pointIndexHit>& info,
1859     labelList& region
1860 ) const
1861 {
1862     if (!Pstream::parRun())
1863     {
1864         region.setSize(info.size());
1865         forAll(info, i)
1866         {
1867             if (info[i].hit())
1868             {
1869                 region[i] = triSurface::operator[](info[i].index()).region();
1870             }
1871             else
1872             {
1873                 region[i] = -1;
1874             }
1875         }
1876         return;
1877     }
1878
1879
1880     // Get query data (= local index of triangle)
1881     // ~~~~~~~~~~~~~~
1882
1883     labelList triangleIndex(info.size());
1884     autoPtr<mapDistribute> mapPtr
1885     (
1886         calcLocalQueries
1887         (
1888             info,
1889             triangleIndex
1890         )
1891     );
1892     const mapDistribute& map = mapPtr();
1893
1894
1895     // Do my tests
1896     // ~~~~~~~~~~~
1897
1898     const triSurface& s = static_cast<const triSurface&>(*this);
1899
1900     region.setSize(triangleIndex.size());
1901
1902     forAll(triangleIndex, i)
1903     {
1904         label triI = triangleIndex[i];
1905         region[i] = s[triI].region();
1906     }
1907
1908
1909     // Send back results
1910     // ~~~~~~~~~~~~~~~~~
1911
1912     map.distribute
1913     (
1914         //Pstream::scheduled,
1915         //map.schedule            // note reverse schedule
1916         //(
1917         //    map.constructMap(),
1918         //    map.subMap()
1919         //),
1920         Pstream::nonBlocking,
1921         List<labelPair>(0),
1922         info.size(),
1923         map.constructMap(),     // what to send
1924         map.subMap(),           // what to receive
1925         region
1926     );
1927 }
1928
1929
1930 void Foam::distributedTriSurfaceMesh::getNormal
1931 (
1932     const List<pointIndexHit>& info,
1933     vectorField& normal
1934 ) const
1935 {
1936     if (!Pstream::parRun())
1937     {
1938         normal.setSize(info.size());
1939
1940         forAll(info, i)
1941         {
1942             if (info[i].hit())
1943             {
1944                 normal[i] = faceNormals()[info[i].index()];
1945             }
1946             else
1947             {
1948                 // Set to what?
1949                 normal[i] = vector::zero;
1950             }
1951         }
1952         return;
1953     }
1954
1955
1956     // Get query data (= local index of triangle)
1957     // ~~~~~~~~~~~~~~
1958
1959     labelList triangleIndex(info.size());
1960     autoPtr<mapDistribute> mapPtr
1961     (
1962         calcLocalQueries
1963         (
1964             info,
1965             triangleIndex
1966         )
1967     );
1968     const mapDistribute& map = mapPtr();
1969
1970
1971     // Do my tests
1972     // ~~~~~~~~~~~
1973
1974     const triSurface& s = static_cast<const triSurface&>(*this);
1975
1976     normal.setSize(triangleIndex.size());
1977
1978     forAll(triangleIndex, i)
1979     {
1980         label triI = triangleIndex[i];
1981         normal[i] = s[triI].normal(s.points());
1982         normal[i] /= mag(normal[i]) + VSMALL;
1983     }
1984
1985
1986     // Send back results
1987     // ~~~~~~~~~~~~~~~~~
1988
1989     map.distribute
1990     (
1991         //Pstream::scheduled,
1992         //map.schedule            // note reverse schedule
1993         //(
1994         //    map.constructMap(),
1995         //    map.subMap()
1996         //),
1997         Pstream::nonBlocking,
1998         List<labelPair>(0),
1999         info.size(),
2000         map.constructMap(),     // what to send
2001         map.subMap(),           // what to receive
2002         normal
2003     );
2004 }
2005
2006
2007 void Foam::distributedTriSurfaceMesh::getField
2008 (
2009     const word& fieldName,
2010     const List<pointIndexHit>& info,
2011     labelList& values
2012 ) const
2013 {
2014     const triSurfaceLabelField& fld = lookupObject<triSurfaceLabelField>
2015     (
2016         fieldName
2017     );
2018
2019
2020     if (!Pstream::parRun())
2021     {
2022         values.setSize(info.size());
2023         forAll(info, i)
2024         {
2025             if (info[i].hit())
2026             {
2027                 values[i] = fld[info[i].index()];
2028             }
2029         }
2030         return;
2031     }
2032
2033
2034     // Get query data (= local index of triangle)
2035     // ~~~~~~~~~~~~~~
2036
2037     labelList triangleIndex(info.size());
2038     autoPtr<mapDistribute> mapPtr
2039     (
2040         calcLocalQueries
2041         (
2042             info,
2043             triangleIndex
2044         )
2045     );
2046     const mapDistribute& map = mapPtr();
2047
2048
2049     // Do my tests
2050     // ~~~~~~~~~~~
2051
2052     values.setSize(triangleIndex.size());
2053
2054     forAll(triangleIndex, i)
2055     {
2056         label triI = triangleIndex[i];
2057         values[i] = fld[triI];
2058     }
2059
2060
2061     // Send back results
2062     // ~~~~~~~~~~~~~~~~~
2063
2064     map.distribute
2065     (
2066         Pstream::nonBlocking,
2067         List<labelPair>(0),
2068         info.size(),
2069         map.constructMap(),     // what to send
2070         map.subMap(),           // what to receive
2071         values
2072     );
2073 }
2074
2075
2076 void Foam::distributedTriSurfaceMesh::getVolumeType
2077 (
2078     const pointField& points,
2079     List<volumeType>& volType
2080 ) const
2081 {
2082     FatalErrorIn
2083     (
2084         "distributedTriSurfaceMesh::getVolumeType"
2085         "(const pointField&, List<volumeType>&) const"
2086     )   << "Volume type not supported for distributed surfaces."
2087         << exit(FatalError);
2088 }
2089
2090
2091 // Subset the part of surface that is overlapping bb.
2092 Foam::triSurface Foam::distributedTriSurfaceMesh::overlappingSurface
2093 (
2094     const triSurface& s,
2095     const List<treeBoundBox>& bbs,
2096
2097     labelList& subPointMap,
2098     labelList& subFaceMap
2099 )
2100 {
2101     // Determine what triangles to keep.
2102     boolList includedFace(s.size(), false);
2103
2104     // Create a slightly larger bounding box.
2105     List<treeBoundBox> bbsX(bbs.size());
2106     const scalar eps = 1.0e-4;
2107     forAll(bbs, i)
2108     {
2109         const point mid = 0.5*(bbs[i].min() + bbs[i].max());
2110         const vector halfSpan = (1.0+eps)*(bbs[i].max() - mid);
2111
2112         bbsX[i].min() = mid - halfSpan;
2113         bbsX[i].max() = mid + halfSpan;
2114     }
2115
2116     forAll(s, triI)
2117     {
2118         const labelledTri& f = s[triI];
2119         const point& p0 = s.points()[f[0]];
2120         const point& p1 = s.points()[f[1]];
2121         const point& p2 = s.points()[f[2]];
2122
2123         if (overlaps(bbsX, p0, p1, p2))
2124         {
2125             includedFace[triI] = true;
2126         }
2127     }
2128
2129     return subsetMesh(s, includedFace, subPointMap, subFaceMap);
2130 }
2131
2132
2133 void Foam::distributedTriSurfaceMesh::distribute
2134 (
2135     const List<treeBoundBox>& bbs,
2136     const bool keepNonLocal,
2137     autoPtr<mapDistribute>& faceMap,
2138     autoPtr<mapDistribute>& pointMap
2139 )
2140 {
2141     // Get bbs of all domains
2142     // ~~~~~~~~~~~~~~~~~~~~~~
2143
2144     {
2145         List<List<treeBoundBox> > newProcBb(Pstream::nProcs());
2146
2147         switch(distType_)
2148         {
2149             case FOLLOW:
2150                 newProcBb[Pstream::myProcNo()].setSize(bbs.size());
2151                 forAll(bbs, i)
2152                 {
2153                     newProcBb[Pstream::myProcNo()][i] = bbs[i];
2154                 }
2155                 Pstream::gatherList(newProcBb);
2156                 Pstream::scatterList(newProcBb);
2157             break;
2158
2159             case INDEPENDENT:
2160                 newProcBb = independentlyDistributedBbs(*this);
2161             break;
2162
2163             case FROZEN:
2164                 return;
2165             break;
2166
2167             default:
2168                 FatalErrorIn("distributedTriSurfaceMesh::distribute(..)")
2169                     << "Unsupported distribution type." << exit(FatalError);
2170             break;
2171         }
2172
2173         //if (debug)
2174         //{
2175         //    Info<< "old bb:" << procBb_ << endl << endl;
2176         //    Info<< "new bb:" << newProcBb << endl << endl;
2177         //    Info<< "Same:" << (newProcBb == procBb_) << endl;
2178         //}
2179
2180         if (newProcBb == procBb_)
2181         {
2182             return;
2183         }
2184         else
2185         {
2186             procBb_.transfer(newProcBb);
2187             dict_.set("bounds", procBb_[Pstream::myProcNo()]);
2188         }
2189     }
2190
2191
2192     // Debug information
2193     if (debug)
2194     {
2195         labelList nTris(Pstream::nProcs());
2196         nTris[Pstream::myProcNo()] = triSurface::size();
2197         Pstream::gatherList(nTris);
2198         Pstream::scatterList(nTris);
2199
2200         Info<< "distributedTriSurfaceMesh::distribute : before distribution:"
2201             << endl
2202             << "\tproc\ttris" << endl;
2203
2204         forAll(nTris, procI)
2205         {
2206             Info<< '\t' << procI << '\t' << nTris[procI] << endl;
2207         }
2208         Info<< endl;
2209     }
2210
2211
2212     // Use procBbs to determine which faces go where
2213     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2214
2215     labelListList faceSendMap(Pstream::nProcs());
2216     labelListList pointSendMap(Pstream::nProcs());
2217
2218     forAll(procBb_, procI)
2219     {
2220         overlappingSurface
2221         (
2222             *this,
2223             procBb_[procI],
2224             pointSendMap[procI],
2225             faceSendMap[procI]
2226         );
2227
2228         if (debug)
2229         {
2230             //Pout<< "Overlapping with proc " << procI
2231             //    << " faces:" << faceSendMap[procI].size()
2232             //    << " points:" << pointSendMap[procI].size() << endl << endl;
2233         }
2234     }
2235
2236     if (keepNonLocal)
2237     {
2238         // Include in faceSendMap/pointSendMap the triangles that are
2239         // not mapped to any processor so they stay local.
2240
2241         const triSurface& s = static_cast<const triSurface&>(*this);
2242
2243         boolList includedFace(s.size(), true);
2244
2245         forAll(faceSendMap, procI)
2246         {
2247             if (procI != Pstream::myProcNo())
2248             {
2249                 forAll(faceSendMap[procI], i)
2250                 {
2251                     includedFace[faceSendMap[procI][i]] = false;
2252                 }
2253             }
2254         }
2255
2256         // Combine includedFace (all triangles that are not on any neighbour)
2257         // with overlap.
2258
2259         forAll(faceSendMap[Pstream::myProcNo()], i)
2260         {
2261             includedFace[faceSendMap[Pstream::myProcNo()][i]] = true;
2262         }
2263
2264         subsetMesh
2265         (
2266             s,
2267             includedFace,
2268             pointSendMap[Pstream::myProcNo()],
2269             faceSendMap[Pstream::myProcNo()]
2270         );
2271     }
2272
2273
2274     // Send over how many faces/points I need to receive
2275     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2276
2277     labelListList faceSendSizes(Pstream::nProcs());
2278     faceSendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
2279     forAll(faceSendMap, procI)
2280     {
2281         faceSendSizes[Pstream::myProcNo()][procI] = faceSendMap[procI].size();
2282     }
2283     Pstream::gatherList(faceSendSizes);
2284     Pstream::scatterList(faceSendSizes);
2285
2286
2287
2288     // Exchange surfaces
2289     // ~~~~~~~~~~~~~~~~~
2290
2291     // Storage for resulting surface
2292     List<labelledTri> allTris;
2293     pointField allPoints;
2294
2295     labelListList faceConstructMap(Pstream::nProcs());
2296     labelListList pointConstructMap(Pstream::nProcs());
2297
2298
2299     // My own surface first
2300     // ~~~~~~~~~~~~~~~~~~~~
2301
2302     {
2303         labelList pointMap;
2304         triSurface subSurface
2305         (
2306             subsetMesh
2307             (
2308                 *this,
2309                 faceSendMap[Pstream::myProcNo()],
2310                 pointMap
2311             )
2312         );
2313
2314         allTris = subSurface;
2315         allPoints = subSurface.points();
2316
2317         faceConstructMap[Pstream::myProcNo()] = identity
2318         (
2319             faceSendMap[Pstream::myProcNo()].size()
2320         );
2321         pointConstructMap[Pstream::myProcNo()] = identity
2322         (
2323             pointSendMap[Pstream::myProcNo()].size()
2324         );
2325     }
2326
2327
2328
2329     // Send all
2330     // ~~~~~~~~
2331
2332     forAll(faceSendSizes, procI)
2333     {
2334         if (procI != Pstream::myProcNo())
2335         {
2336             if (faceSendSizes[Pstream::myProcNo()][procI] > 0)
2337             {
2338                 OPstream str(Pstream::blocking, procI);
2339
2340                 labelList pointMap;
2341                 triSurface subSurface
2342                 (
2343                     subsetMesh
2344                     (
2345                         *this,
2346                         faceSendMap[procI],
2347                         pointMap
2348                     )
2349                 );
2350
2351                 //if (debug)
2352                 //{
2353                 //    Pout<< "Sending to " << procI
2354                 //        << " faces:" << faceSendMap[procI].size()
2355                 //        << " points:" << subSurface.points().size() << endl
2356                 //        << endl;
2357                 //}
2358
2359                 str << subSurface;
2360            }
2361         }
2362     }
2363
2364
2365     // Receive and merge all
2366     // ~~~~~~~~~~~~~~~~~~~~~
2367
2368     forAll(faceSendSizes, procI)
2369     {
2370         if (procI != Pstream::myProcNo())
2371         {
2372             if (faceSendSizes[procI][Pstream::myProcNo()] > 0)
2373             {
2374                 IPstream str(Pstream::blocking, procI);
2375
2376                 // Receive
2377                 triSurface subSurface(str);
2378
2379                 //if (debug)
2380                 //{
2381                 //    Pout<< "Received from " << procI
2382                 //        << " faces:" << subSurface.size()
2383                 //        << " points:" << subSurface.points().size() << endl
2384                 //        << endl;
2385                 //}
2386
2387                 // Merge into allSurf
2388                 merge
2389                 (
2390                     mergeDist_,
2391                     subSurface,
2392                     subSurface.points(),
2393
2394                     allTris,
2395                     allPoints,
2396                     faceConstructMap[procI],
2397                     pointConstructMap[procI]
2398                 );
2399
2400                 //if (debug)
2401                 //{
2402                 //    Pout<< "Current merged surface : faces:" << allTris.size()
2403                 //        << " points:" << allPoints.size() << endl << endl;
2404                 //}
2405            }
2406         }
2407     }
2408
2409
2410     faceMap.reset
2411     (
2412         new mapDistribute
2413         (
2414             allTris.size(),
2415             faceSendMap,
2416             faceConstructMap,
2417             true
2418         )
2419     );
2420     pointMap.reset
2421     (
2422         new mapDistribute
2423         (
2424             allPoints.size(),
2425             pointSendMap,
2426             pointConstructMap,
2427             true
2428         )
2429     );
2430
2431     // Construct triSurface. Reuse storage.
2432     triSurface::operator=(triSurface(allTris, patches(), allPoints, true));
2433
2434     clearOut();
2435
2436     // Regions stays same
2437     // Volume type stays same.
2438
2439     distributeFields<label>(faceMap());
2440     distributeFields<scalar>(faceMap());
2441     distributeFields<vector>(faceMap());
2442     distributeFields<sphericalTensor>(faceMap());
2443     distributeFields<symmTensor>(faceMap());
2444     distributeFields<tensor>(faceMap());
2445
2446     if (debug)
2447     {
2448         labelList nTris(Pstream::nProcs());
2449         nTris[Pstream::myProcNo()] = triSurface::size();
2450         Pstream::gatherList(nTris);
2451         Pstream::scatterList(nTris);
2452
2453         Info<< "distributedTriSurfaceMesh::distribute : after distribution:"
2454             << endl
2455             << "\tproc\ttris" << endl;
2456
2457         forAll(nTris, procI)
2458         {
2459             Info<< '\t' << procI << '\t' << nTris[procI] << endl;
2460         }
2461         Info<< endl;
2462     }
2463 }
2464
2465
2466 //- Write using given format, version and compression
2467 bool Foam::distributedTriSurfaceMesh::writeObject
2468 (
2469     IOstream::streamFormat fmt,
2470     IOstream::versionNumber ver,
2471     IOstream::compressionType cmp
2472 ) const
2473 {
2474     // Make sure dictionary goes to same directory as surface
2475     const_cast<fileName&>(dict_.instance()) = searchableSurface::instance();
2476
2477     // Dictionary needs to be written in ascii - binary output not supported.
2478     return
2479         triSurfaceMesh::writeObject(fmt, ver, cmp)
2480      && dict_.writeObject(IOstream::ASCII, ver, cmp);
2481 }
2482
2483
2484 void Foam::distributedTriSurfaceMesh::writeStats(Ostream& os) const
2485 {
2486     boundBox bb;
2487     label nPoints;
2488     calcBounds(bb, nPoints);
2489     reduce(bb.min(), minOp<point>());
2490     reduce(bb.max(), maxOp<point>());
2491
2492     os  << "Triangles    : " << returnReduce(triSurface::size(), sumOp<label>())
2493         << endl
2494         << "Vertices     : " << returnReduce(nPoints, sumOp<label>()) << endl
2495         << "Bounding Box : " << bb << endl;
2496 }
2497
2498
2499 // ************************************************************************* //