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