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_(io, dict)
1340 {
1341     read();
1342
1343     if (debug)
1344     {
1345         Info<< "Constructed from triSurface:" << endl;
1346         writeStats(Info);
1347
1348         labelList nTris(Pstream::nProcs());
1349         nTris[Pstream::myProcNo()] = triSurface::size();
1350         Pstream::gatherList(nTris);
1351         Pstream::scatterList(nTris);
1352
1353         Info<< endl<< "\tproc\ttris\tbb" << endl;
1354         forAll(nTris, procI)
1355         {
1356             Info<< '\t' << procI << '\t' << nTris[procI]
1357                  << '\t' << procBb_[procI] << endl;
1358         }
1359         Info<< endl;
1360     }
1361 }
1362
1363
1364 Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh(const IOobject& io)
1365 :
1366     //triSurfaceMesh(io),
1367     triSurfaceMesh
1368     (
1369         IOobject
1370         (
1371             io.name(),
1372             io.time().findInstance(io.local(), word::null),
1373             io.local(),
1374             io.db(),
1375             io.readOpt(),
1376             io.writeOpt(),
1377             io.registerObject()
1378         )
1379     ),
1380     dict_
1381     (
1382         IOobject
1383         (
1384             searchableSurface::name() + "Dict",
1385             searchableSurface::instance(),
1386             searchableSurface::local(),
1387             searchableSurface::db(),
1388             searchableSurface::readOpt(),
1389             searchableSurface::writeOpt(),
1390             searchableSurface::registerObject()
1391         )
1392     )
1393 {
1394     read();
1395
1396     if (debug)
1397     {
1398         Info<< "Read distributedTriSurface from " << io.objectPath()
1399             << ':' << endl;
1400         writeStats(Info);
1401
1402         labelList nTris(Pstream::nProcs());
1403         nTris[Pstream::myProcNo()] = triSurface::size();
1404         Pstream::gatherList(nTris);
1405         Pstream::scatterList(nTris);
1406
1407         Info<< endl<< "\tproc\ttris\tbb" << endl;
1408         forAll(nTris, procI)
1409         {
1410             Info<< '\t' << procI << '\t' << nTris[procI]
1411                  << '\t' << procBb_[procI] << endl;
1412         }
1413         Info<< endl;
1414     }
1415 }
1416
1417
1418 Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh
1419 (
1420     const IOobject& io,
1421     const dictionary& dict
1422 )
1423 :
1424     //triSurfaceMesh(io, dict),
1425     triSurfaceMesh
1426     (
1427         IOobject
1428         (
1429             io.name(),
1430             io.time().findInstance(io.local(), word::null),
1431             io.local(),
1432             io.db(),
1433             io.readOpt(),
1434             io.writeOpt(),
1435             io.registerObject()
1436         ),
1437         dict
1438     ),
1439     dict_
1440     (
1441         IOobject
1442         (
1443             searchableSurface::name() + "Dict",
1444             searchableSurface::instance(),
1445             searchableSurface::local(),
1446             searchableSurface::db(),
1447             searchableSurface::readOpt(),
1448             searchableSurface::writeOpt(),
1449             searchableSurface::registerObject()
1450         )
1451     )
1452 {
1453     read();
1454
1455     if (debug)
1456     {
1457         Info<< "Read distributedTriSurface from " << io.objectPath()
1458             << " and dictionary:" << endl;
1459         writeStats(Info);
1460
1461         labelList nTris(Pstream::nProcs());
1462         nTris[Pstream::myProcNo()] = triSurface::size();
1463         Pstream::gatherList(nTris);
1464         Pstream::scatterList(nTris);
1465
1466         Info<< endl<< "\tproc\ttris\tbb" << endl;
1467         forAll(nTris, procI)
1468         {
1469             Info<< '\t' << procI << '\t' << nTris[procI]
1470                  << '\t' << procBb_[procI] << endl;
1471         }
1472         Info<< endl;
1473     }
1474 }
1475
1476
1477 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
1478
1479 Foam::distributedTriSurfaceMesh::~distributedTriSurfaceMesh()
1480 {
1481     clearOut();
1482 }
1483
1484
1485 void Foam::distributedTriSurfaceMesh::clearOut()
1486 {
1487     globalTris_.clear();
1488     triSurfaceMesh::clearOut();
1489 }
1490
1491
1492 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
1493
1494 const Foam::globalIndex& Foam::distributedTriSurfaceMesh::globalTris() const
1495 {
1496     if (!globalTris_.valid())
1497     {
1498         globalTris_.reset(new globalIndex(triSurface::size()));
1499     }
1500     return globalTris_;
1501 }
1502
1503
1504 void Foam::distributedTriSurfaceMesh::findNearest
1505 (
1506     const pointField& samples,
1507     const scalarField& nearestDistSqr,
1508     List<pointIndexHit>& info
1509 ) const
1510 {
1511     const indexedOctree<treeDataTriSurface>& octree = tree();
1512
1513     // Important:force synchronised construction of indexing
1514     const globalIndex& triIndexer = globalTris();
1515
1516
1517     // Initialise
1518     // ~~~~~~~~~~
1519
1520     info.setSize(samples.size());
1521     forAll(info, i)
1522     {
1523         info[i].setMiss();
1524     }
1525
1526
1527
1528     // Do any local queries
1529     // ~~~~~~~~~~~~~~~~~~~~
1530
1531     label nLocal = 0;
1532
1533     {
1534         // Work array - whether processor bb overlaps the bounding sphere.
1535         boolList procBbOverlaps(Pstream::nProcs());
1536
1537         forAll(samples, i)
1538         {
1539             // Find the processor this sample+radius overlaps.
1540             label nProcs = calcOverlappingProcs
1541             (
1542                 samples[i],
1543                 nearestDistSqr[i],
1544                 procBbOverlaps
1545             );
1546
1547             // Overlaps local processor?
1548             if (procBbOverlaps[Pstream::myProcNo()])
1549             {
1550                 info[i] = octree.findNearest(samples[i], nearestDistSqr[i]);
1551                 if (info[i].hit())
1552                 {
1553                     info[i].setIndex(triIndexer.toGlobal(info[i].index()));
1554                 }
1555                 if (nProcs == 1)
1556                 {
1557                     // Fully local
1558                     nLocal++;
1559                 }
1560             }
1561         }
1562     }
1563
1564
1565     if
1566     (
1567         Pstream::parRun()
1568      && (
1569             returnReduce(nLocal, sumOp<label>())
1570           < returnReduce(samples.size(), sumOp<label>())
1571         )
1572     )
1573     {
1574         // Not all can be resolved locally. Build queries and map, send over
1575         // queries, do intersections, send back and merge.
1576
1577         // Calculate queries and exchange map
1578         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1579
1580         pointField allCentres;
1581         scalarField allRadiusSqr;
1582         labelList allSegmentMap;
1583         autoPtr<mapDistribute> mapPtr
1584         (
1585             calcLocalQueries
1586             (
1587                 samples,
1588                 nearestDistSqr,
1589
1590                 allCentres,
1591                 allRadiusSqr,
1592                 allSegmentMap
1593             )
1594         );
1595         const mapDistribute& map = mapPtr();
1596
1597
1598         // swap samples to local processor
1599         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1600
1601         map.distribute
1602         (
1603             //Pstream::scheduled,
1604             //map.schedule(),
1605             Pstream::nonBlocking,
1606             List<labelPair>(0),
1607             map.constructSize(),
1608             map.subMap(),           // what to send
1609             map.constructMap(),     // what to receive
1610             allCentres
1611         );
1612         map.distribute
1613         (
1614             //Pstream::scheduled,
1615             //map.schedule(),
1616             Pstream::nonBlocking,
1617             List<labelPair>(0),
1618             map.constructSize(),
1619             map.subMap(),           // what to send
1620             map.constructMap(),     // what to receive
1621             allRadiusSqr
1622         );
1623
1624
1625         // Do my tests
1626         // ~~~~~~~~~~~
1627
1628         List<pointIndexHit> allInfo(allCentres.size());
1629         forAll(allInfo, i)
1630         {
1631             allInfo[i] = octree.findNearest
1632             (
1633                 allCentres[i],
1634                 allRadiusSqr[i]
1635             );
1636             if (allInfo[i].hit())
1637             {
1638                 allInfo[i].setIndex(triIndexer.toGlobal(allInfo[i].index()));
1639             }
1640         }
1641
1642
1643         // Send back results
1644         // ~~~~~~~~~~~~~~~~~
1645
1646         map.distribute
1647         (
1648             //Pstream::scheduled,
1649             //map.schedule            // note reverse schedule
1650             //(
1651             //    map.constructMap(),
1652             //    map.subMap()
1653             //),
1654             Pstream::nonBlocking,
1655             List<labelPair>(0),
1656             allSegmentMap.size(),
1657             map.constructMap(),     // what to send
1658             map.subMap(),           // what to receive
1659             allInfo
1660         );
1661
1662
1663         // Extract information
1664         // ~~~~~~~~~~~~~~~~~~~
1665
1666         forAll(allInfo, i)
1667         {
1668             if (allInfo[i].hit())
1669             {
1670                 label pointI = allSegmentMap[i];
1671
1672                 if (!info[pointI].hit())
1673                 {
1674                     // No intersection yet so take this one
1675                     info[pointI] = allInfo[i];
1676                 }
1677                 else
1678                 {
1679                     // Nearest intersection
1680                     if
1681                     (
1682                         magSqr(allInfo[i].hitPoint()-samples[pointI])
1683                       < magSqr(info[pointI].hitPoint()-samples[pointI])
1684                     )
1685                     {
1686                         info[pointI] = allInfo[i];
1687                     }
1688                 }
1689             }
1690         }
1691     }
1692 }
1693
1694
1695 void Foam::distributedTriSurfaceMesh::findLine
1696 (
1697     const pointField& start,
1698     const pointField& end,
1699     List<pointIndexHit>& info
1700 ) const
1701 {
1702     findLine
1703     (
1704         true,   // nearestIntersection
1705         start,
1706         end,
1707         info
1708     );
1709 }
1710
1711
1712 void Foam::distributedTriSurfaceMesh::findLineAny
1713 (
1714     const pointField& start,
1715     const pointField& end,
1716     List<pointIndexHit>& info
1717 ) const
1718 {
1719     findLine
1720     (
1721         true,   // nearestIntersection
1722         start,
1723         end,
1724         info
1725     );
1726 }
1727
1728
1729 void Foam::distributedTriSurfaceMesh::findLineAll
1730 (
1731     const pointField& start,
1732     const pointField& end,
1733     List<List<pointIndexHit> >& info
1734 ) const
1735 {
1736     // Reuse fineLine. We could modify all of findLine to do multiple
1737     // intersections but this would mean a lot of data transferred so
1738     // for now we just find nearest intersection and retest from that
1739     // intersection onwards.
1740
1741     // Work array.
1742     List<pointIndexHit> hitInfo(start.size());
1743
1744     findLine
1745     (
1746         true,   // nearestIntersection
1747         start,
1748         end,
1749         hitInfo
1750     );
1751
1752     // Tolerances:
1753     // To find all intersections we add a small vector to the last intersection
1754     // This is chosen such that
1755     // - it is significant (SMALL is smallest representative relative tolerance;
1756     //   we need something bigger since we're doing calculations)
1757     // - if the start-end vector is zero we still progress
1758     const vectorField dirVec(end-start);
1759     const scalarField magSqrDirVec(magSqr(dirVec));
1760     const vectorField smallVec
1761     (
1762         Foam::sqrt(SMALL)*dirVec
1763       + vector(ROOTVSMALL,ROOTVSMALL,ROOTVSMALL)
1764     );
1765
1766     // Copy to input and compact any hits
1767     labelList pointMap(start.size());
1768     pointField e0(start.size());
1769     pointField e1(start.size());
1770     label compactI = 0;
1771
1772     info.setSize(hitInfo.size());
1773     forAll(hitInfo, pointI)
1774     {
1775         if (hitInfo[pointI].hit())
1776         {
1777             info[pointI].setSize(1);
1778             info[pointI][0] = hitInfo[pointI];
1779
1780             point pt = hitInfo[pointI].hitPoint() + smallVec[pointI];
1781
1782             if (((pt-start[pointI])&dirVec[pointI]) <= magSqrDirVec[pointI])
1783             {
1784                 e0[compactI] = pt;
1785                 e1[compactI] = end[pointI];
1786                 pointMap[compactI] = pointI;
1787                 compactI++;
1788             }
1789         }
1790         else
1791         {
1792             info[pointI].clear();
1793         }
1794     }
1795
1796     e0.setSize(compactI);
1797     e1.setSize(compactI);
1798     pointMap.setSize(compactI);
1799
1800     while (returnReduce(e0.size(), sumOp<label>()) > 0)
1801     {
1802         findLine
1803         (
1804             true,   // nearestIntersection
1805             e0,
1806             e1,
1807             hitInfo
1808         );
1809
1810
1811         // Extract
1812         label compactI = 0;
1813         forAll(hitInfo, i)
1814         {
1815             if (hitInfo[i].hit())
1816             {
1817                 label pointI = pointMap[i];
1818
1819                 label sz = info[pointI].size();
1820                 info[pointI].setSize(sz+1);
1821                 info[pointI][sz] = hitInfo[i];
1822
1823                 point pt = hitInfo[i].hitPoint() + smallVec[pointI];
1824
1825                 if (((pt-start[pointI])&dirVec[pointI]) <= magSqrDirVec[pointI])
1826                 {
1827                     e0[compactI] = pt;
1828                     e1[compactI] = end[pointI];
1829                     pointMap[compactI] = pointI;
1830                     compactI++;
1831                 }
1832             }
1833         }
1834
1835         // Trim
1836         e0.setSize(compactI);
1837         e1.setSize(compactI);
1838         pointMap.setSize(compactI);
1839     }
1840 }
1841
1842
1843 void Foam::distributedTriSurfaceMesh::getRegion
1844 (
1845     const List<pointIndexHit>& info,
1846     labelList& region
1847 ) const
1848 {
1849     if (!Pstream::parRun())
1850     {
1851         region.setSize(info.size());
1852         forAll(info, i)
1853         {
1854             if (info[i].hit())
1855             {
1856                 region[i] = triSurface::operator[](info[i].index()).region();
1857             }
1858             else
1859             {
1860                 region[i] = -1;
1861             }
1862         }
1863         return;
1864     }
1865
1866
1867     // Get query data (= local index of triangle)
1868     // ~~~~~~~~~~~~~~
1869
1870     labelList triangleIndex(info.size());
1871     autoPtr<mapDistribute> mapPtr
1872     (
1873         calcLocalQueries
1874         (
1875             info,
1876             triangleIndex
1877         )
1878     );
1879     const mapDistribute& map = mapPtr();
1880
1881
1882     // Do my tests
1883     // ~~~~~~~~~~~
1884
1885     const triSurface& s = static_cast<const triSurface&>(*this);
1886
1887     region.setSize(triangleIndex.size());
1888
1889     forAll(triangleIndex, i)
1890     {
1891         label triI = triangleIndex[i];
1892         region[i] = s[triI].region();
1893     }
1894
1895
1896     // Send back results
1897     // ~~~~~~~~~~~~~~~~~
1898
1899     map.distribute
1900     (
1901         //Pstream::scheduled,
1902         //map.schedule            // note reverse schedule
1903         //(
1904         //    map.constructMap(),
1905         //    map.subMap()
1906         //),
1907         Pstream::nonBlocking,
1908         List<labelPair>(0),
1909         info.size(),
1910         map.constructMap(),     // what to send
1911         map.subMap(),           // what to receive
1912         region
1913     );
1914 }
1915
1916
1917 void Foam::distributedTriSurfaceMesh::getNormal
1918 (
1919     const List<pointIndexHit>& info,
1920     vectorField& normal
1921 ) const
1922 {
1923     if (!Pstream::parRun())
1924     {
1925         normal.setSize(info.size());
1926
1927         forAll(info, i)
1928         {
1929             if (info[i].hit())
1930             {
1931                 normal[i] = faceNormals()[info[i].index()];
1932             }
1933             else
1934             {
1935                 // Set to what?
1936                 normal[i] = vector::zero;
1937             }
1938         }
1939         return;
1940     }
1941
1942
1943     // Get query data (= local index of triangle)
1944     // ~~~~~~~~~~~~~~
1945
1946     labelList triangleIndex(info.size());
1947     autoPtr<mapDistribute> mapPtr
1948     (
1949         calcLocalQueries
1950         (
1951             info,
1952             triangleIndex
1953         )
1954     );
1955     const mapDistribute& map = mapPtr();
1956
1957
1958     // Do my tests
1959     // ~~~~~~~~~~~
1960
1961     const triSurface& s = static_cast<const triSurface&>(*this);
1962
1963     normal.setSize(triangleIndex.size());
1964
1965     forAll(triangleIndex, i)
1966     {
1967         label triI = triangleIndex[i];
1968         normal[i] = s[triI].normal(s.points());
1969         normal[i] /= mag(normal[i]) + VSMALL;
1970     }
1971
1972
1973     // Send back results
1974     // ~~~~~~~~~~~~~~~~~
1975
1976     map.distribute
1977     (
1978         //Pstream::scheduled,
1979         //map.schedule            // note reverse schedule
1980         //(
1981         //    map.constructMap(),
1982         //    map.subMap()
1983         //),
1984         Pstream::nonBlocking,
1985         List<labelPair>(0),
1986         info.size(),
1987         map.constructMap(),     // what to send
1988         map.subMap(),           // what to receive
1989         normal
1990     );
1991 }
1992
1993
1994 void Foam::distributedTriSurfaceMesh::getField
1995 (
1996     const word& fieldName,
1997     const List<pointIndexHit>& info,
1998     labelList& values
1999 ) const
2000 {
2001     const triSurfaceLabelField& fld = lookupObject<triSurfaceLabelField>
2002     (
2003         fieldName
2004     );
2005
2006
2007     if (!Pstream::parRun())
2008     {
2009         values.setSize(info.size());
2010         forAll(info, i)
2011         {
2012             if (info[i].hit())
2013             {
2014                 values[i] = fld[info[i].index()];
2015             }
2016         }
2017         return;
2018     }
2019
2020
2021     // Get query data (= local index of triangle)
2022     // ~~~~~~~~~~~~~~
2023
2024     labelList triangleIndex(info.size());
2025     autoPtr<mapDistribute> mapPtr
2026     (
2027         calcLocalQueries
2028         (
2029             info,
2030             triangleIndex
2031         )
2032     );
2033     const mapDistribute& map = mapPtr();
2034
2035
2036     // Do my tests
2037     // ~~~~~~~~~~~
2038
2039     values.setSize(triangleIndex.size());
2040
2041     forAll(triangleIndex, i)
2042     {
2043         label triI = triangleIndex[i];
2044         values[i] = fld[triI];
2045     }
2046
2047
2048     // Send back results
2049     // ~~~~~~~~~~~~~~~~~
2050
2051     map.distribute
2052     (
2053         Pstream::nonBlocking,
2054         List<labelPair>(0),
2055         info.size(),
2056         map.constructMap(),     // what to send
2057         map.subMap(),           // what to receive
2058         values
2059     );
2060 }
2061
2062
2063 void Foam::distributedTriSurfaceMesh::getVolumeType
2064 (
2065     const pointField& points,
2066     List<volumeType>& volType
2067 ) const
2068 {
2069     FatalErrorIn
2070     (
2071         "distributedTriSurfaceMesh::getVolumeType"
2072         "(const pointField&, List<volumeType>&) const"
2073     )   << "Volume type not supported for distributed surfaces."
2074         << exit(FatalError);
2075 }
2076
2077
2078 // Subset the part of surface that is overlapping bb.
2079 Foam::triSurface Foam::distributedTriSurfaceMesh::overlappingSurface
2080 (
2081     const triSurface& s,
2082     const List<treeBoundBox>& bbs,
2083
2084     labelList& subPointMap,
2085     labelList& subFaceMap
2086 )
2087 {
2088     // Determine what triangles to keep.
2089     boolList includedFace(s.size(), false);
2090
2091     // Create a slightly larger bounding box.
2092     List<treeBoundBox> bbsX(bbs.size());
2093     const scalar eps = 1.0e-4;
2094     forAll(bbs, i)
2095     {
2096         const point mid = 0.5*(bbs[i].min() + bbs[i].max());
2097         const vector halfSpan = (1.0+eps)*(bbs[i].max() - mid);
2098
2099         bbsX[i].min() = mid - halfSpan;
2100         bbsX[i].max() = mid + halfSpan;
2101     }
2102
2103     forAll(s, triI)
2104     {
2105         const labelledTri& f = s[triI];
2106         const point& p0 = s.points()[f[0]];
2107         const point& p1 = s.points()[f[1]];
2108         const point& p2 = s.points()[f[2]];
2109
2110         if (overlaps(bbsX, p0, p1, p2))
2111         {
2112             includedFace[triI] = true;
2113         }
2114     }
2115
2116     return subsetMesh(s, includedFace, subPointMap, subFaceMap);
2117 }
2118
2119
2120 void Foam::distributedTriSurfaceMesh::distribute
2121 (
2122     const List<treeBoundBox>& bbs,
2123     const bool keepNonLocal,
2124     autoPtr<mapDistribute>& faceMap,
2125     autoPtr<mapDistribute>& pointMap
2126 )
2127 {
2128     // Get bbs of all domains
2129     // ~~~~~~~~~~~~~~~~~~~~~~
2130
2131     {
2132         List<List<treeBoundBox> > newProcBb(Pstream::nProcs());
2133
2134         switch(distType_)
2135         {
2136             case FOLLOW:
2137                 newProcBb[Pstream::myProcNo()].setSize(bbs.size());
2138                 forAll(bbs, i)
2139                 {
2140                     newProcBb[Pstream::myProcNo()][i] = bbs[i];
2141                 }
2142                 Pstream::gatherList(newProcBb);
2143                 Pstream::scatterList(newProcBb);
2144             break;
2145
2146             case INDEPENDENT:
2147                 newProcBb = independentlyDistributedBbs(*this);
2148             break;
2149
2150             case FROZEN:
2151                 return;
2152             break;
2153
2154             default:
2155                 FatalErrorIn("distributedTriSurfaceMesh::distribute(..)")
2156                     << "Unsupported distribution type." << exit(FatalError);
2157             break;
2158         }
2159
2160         //if (debug)
2161         //{
2162         //    Info<< "old bb:" << procBb_ << endl << endl;
2163         //    Info<< "new bb:" << newProcBb << endl << endl;
2164         //    Info<< "Same:" << (newProcBb == procBb_) << endl;
2165         //}
2166
2167         if (newProcBb == procBb_)
2168         {
2169             return;
2170         }
2171         else
2172         {
2173             procBb_.transfer(newProcBb);
2174             dict_.set("bounds", procBb_[Pstream::myProcNo()]);
2175         }
2176     }
2177
2178
2179     // Debug information
2180     if (debug)
2181     {
2182         labelList nTris(Pstream::nProcs());
2183         nTris[Pstream::myProcNo()] = triSurface::size();
2184         Pstream::gatherList(nTris);
2185         Pstream::scatterList(nTris);
2186
2187         Info<< "distributedTriSurfaceMesh::distribute : before distribution:"
2188             << endl
2189             << "\tproc\ttris" << endl;
2190
2191         forAll(nTris, procI)
2192         {
2193             Info<< '\t' << procI << '\t' << nTris[procI] << endl;
2194         }
2195         Info<< endl;
2196     }
2197
2198
2199     // Use procBbs to determine which faces go where
2200     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2201
2202     labelListList faceSendMap(Pstream::nProcs());
2203     labelListList pointSendMap(Pstream::nProcs());
2204
2205     forAll(procBb_, procI)
2206     {
2207         overlappingSurface
2208         (
2209             *this,
2210             procBb_[procI],
2211             pointSendMap[procI],
2212             faceSendMap[procI]
2213         );
2214
2215         if (debug)
2216         {
2217             //Pout<< "Overlapping with proc " << procI
2218             //    << " faces:" << faceSendMap[procI].size()
2219             //    << " points:" << pointSendMap[procI].size() << endl << endl;
2220         }
2221     }
2222
2223     if (keepNonLocal)
2224     {
2225         // Include in faceSendMap/pointSendMap the triangles that are
2226         // not mapped to any processor so they stay local.
2227
2228         const triSurface& s = static_cast<const triSurface&>(*this);
2229
2230         boolList includedFace(s.size(), true);
2231
2232         forAll(faceSendMap, procI)
2233         {
2234             if (procI != Pstream::myProcNo())
2235             {
2236                 forAll(faceSendMap[procI], i)
2237                 {
2238                     includedFace[faceSendMap[procI][i]] = false;
2239                 }
2240             }
2241         }
2242
2243         // Combine includedFace (all triangles that are not on any neighbour)
2244         // with overlap.
2245
2246         forAll(faceSendMap[Pstream::myProcNo()], i)
2247         {
2248             includedFace[faceSendMap[Pstream::myProcNo()][i]] = true;
2249         }
2250
2251         subsetMesh
2252         (
2253             s,
2254             includedFace,
2255             pointSendMap[Pstream::myProcNo()],
2256             faceSendMap[Pstream::myProcNo()]
2257         );
2258     }
2259
2260
2261     // Send over how many faces/points I need to receive
2262     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2263
2264     labelListList faceSendSizes(Pstream::nProcs());
2265     faceSendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
2266     forAll(faceSendMap, procI)
2267     {
2268         faceSendSizes[Pstream::myProcNo()][procI] = faceSendMap[procI].size();
2269     }
2270     Pstream::gatherList(faceSendSizes);
2271     Pstream::scatterList(faceSendSizes);
2272
2273
2274
2275     // Exchange surfaces
2276     // ~~~~~~~~~~~~~~~~~
2277
2278     // Storage for resulting surface
2279     List<labelledTri> allTris;
2280     pointField allPoints;
2281
2282     labelListList faceConstructMap(Pstream::nProcs());
2283     labelListList pointConstructMap(Pstream::nProcs());
2284
2285
2286     // My own surface first
2287     // ~~~~~~~~~~~~~~~~~~~~
2288
2289     {
2290         labelList pointMap;
2291         triSurface subSurface
2292         (
2293             subsetMesh
2294             (
2295                 *this,
2296                 faceSendMap[Pstream::myProcNo()],
2297                 pointMap
2298             )
2299         );
2300
2301         allTris = subSurface;
2302         allPoints = subSurface.points();
2303
2304         faceConstructMap[Pstream::myProcNo()] = identity
2305         (
2306             faceSendMap[Pstream::myProcNo()].size()
2307         );
2308         pointConstructMap[Pstream::myProcNo()] = identity
2309         (
2310             pointSendMap[Pstream::myProcNo()].size()
2311         );
2312     }
2313
2314
2315
2316     // Send all
2317     // ~~~~~~~~
2318
2319     forAll(faceSendSizes, procI)
2320     {
2321         if (procI != Pstream::myProcNo())
2322         {
2323             if (faceSendSizes[Pstream::myProcNo()][procI] > 0)
2324             {
2325                 OPstream str(Pstream::blocking, procI);
2326
2327                 labelList pointMap;
2328                 triSurface subSurface
2329                 (
2330                     subsetMesh
2331                     (
2332                         *this,
2333                         faceSendMap[procI],
2334                         pointMap
2335                     )
2336                 );
2337
2338                 //if (debug)
2339                 //{
2340                 //    Pout<< "Sending to " << procI
2341                 //        << " faces:" << faceSendMap[procI].size()
2342                 //        << " points:" << subSurface.points().size() << endl
2343                 //        << endl;
2344                 //}
2345
2346                 str << subSurface;
2347            }
2348         }
2349     }
2350
2351
2352     // Receive and merge all
2353     // ~~~~~~~~~~~~~~~~~~~~~
2354
2355     forAll(faceSendSizes, procI)
2356     {
2357         if (procI != Pstream::myProcNo())
2358         {
2359             if (faceSendSizes[procI][Pstream::myProcNo()] > 0)
2360             {
2361                 IPstream str(Pstream::blocking, procI);
2362
2363                 // Receive
2364                 triSurface subSurface(str);
2365
2366                 //if (debug)
2367                 //{
2368                 //    Pout<< "Received from " << procI
2369                 //        << " faces:" << subSurface.size()
2370                 //        << " points:" << subSurface.points().size() << endl
2371                 //        << endl;
2372                 //}
2373
2374                 // Merge into allSurf
2375                 merge
2376                 (
2377                     mergeDist_,
2378                     subSurface,
2379                     subSurface.points(),
2380
2381                     allTris,
2382                     allPoints,
2383                     faceConstructMap[procI],
2384                     pointConstructMap[procI]
2385                 );
2386
2387                 //if (debug)
2388                 //{
2389                 //    Pout<< "Current merged surface : faces:" << allTris.size()
2390                 //        << " points:" << allPoints.size() << endl << endl;
2391                 //}
2392            }
2393         }
2394     }
2395
2396
2397     faceMap.reset
2398     (
2399         new mapDistribute
2400         (
2401             allTris.size(),
2402             faceSendMap,
2403             faceConstructMap,
2404             true
2405         )
2406     );
2407     pointMap.reset
2408     (
2409         new mapDistribute
2410         (
2411             allPoints.size(),
2412             pointSendMap,
2413             pointConstructMap,
2414             true
2415         )
2416     );
2417
2418     // Construct triSurface. Reuse storage.
2419     triSurface::operator=(triSurface(allTris, patches(), allPoints, true));
2420
2421     clearOut();
2422
2423     // Regions stays same
2424     // Volume type stays same.
2425
2426     distributeFields<label>(faceMap());
2427     distributeFields<scalar>(faceMap());
2428     distributeFields<vector>(faceMap());
2429     distributeFields<sphericalTensor>(faceMap());
2430     distributeFields<symmTensor>(faceMap());
2431     distributeFields<tensor>(faceMap());
2432
2433     if (debug)
2434     {
2435         labelList nTris(Pstream::nProcs());
2436         nTris[Pstream::myProcNo()] = triSurface::size();
2437         Pstream::gatherList(nTris);
2438         Pstream::scatterList(nTris);
2439
2440         Info<< "distributedTriSurfaceMesh::distribute : after distribution:"
2441             << endl
2442             << "\tproc\ttris" << endl;
2443
2444         forAll(nTris, procI)
2445         {
2446             Info<< '\t' << procI << '\t' << nTris[procI] << endl;
2447         }
2448         Info<< endl;
2449     }
2450 }
2451
2452
2453 //- Write using given format, version and compression
2454 bool Foam::distributedTriSurfaceMesh::writeObject
2455 (
2456     IOstream::streamFormat fmt,
2457     IOstream::versionNumber ver,
2458     IOstream::compressionType cmp
2459 ) const
2460 {
2461     // Make sure dictionary goes to same directory as surface
2462     const_cast<fileName&>(dict_.instance()) = searchableSurface::instance();
2463
2464     // Dictionary needs to be written in ascii - binary output not supported.
2465     return
2466         triSurfaceMesh::writeObject(fmt, ver, cmp)
2467      && dict_.writeObject(IOstream::ASCII, ver, cmp);
2468 }
2469
2470
2471 void Foam::distributedTriSurfaceMesh::writeStats(Ostream& os) const
2472 {
2473     boundBox bb;
2474     label nPoints;
2475     calcBounds(bb, nPoints);
2476     reduce(bb.min(), minOp<point>());
2477     reduce(bb.max(), maxOp<point>());
2478
2479     os  << "Triangles    : " << returnReduce(triSurface::size(), sumOp<label>())
2480         << endl
2481         << "Vertices     : " << returnReduce(nPoints, sumOp<label>()) << endl
2482         << "Bounding Box : " << bb << endl;
2483 }
2484
2485
2486 // ************************************************************************* //