src/meshTools/searchableSurface/distributedTriSurfaceMesh.C

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