FreeFOAM The Cross-Platform CFD Toolkit
Hosted by SourceForge:
Get FreeFOAM at SourceForge.net.
            Fast, secure and Free Open Source software downloads

distributedTriSurfaceMesh.C

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*\
00002   =========                 |
00003   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
00004    \\    /   O peration     |
00005     \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
00006      \\/     M anipulation  |
00007 -------------------------------------------------------------------------------
00008 License
00009     This file is part of OpenFOAM.
00010 
00011     OpenFOAM is free software: you can redistribute it and/or modify it
00012     under the terms of the GNU General Public License as published by
00013     the Free Software Foundation, either version 3 of the License, or
00014     (at your option) any later version.
00015 
00016     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00019     for more details.
00020 
00021     You should have received a copy of the GNU General Public License
00022     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
00023 
00024 \*---------------------------------------------------------------------------*/
00025 
00026 #include "distributedTriSurfaceMesh.H"
00027 #include <OpenFOAM/mapDistribute.H>
00028 #include <OpenFOAM/Random.H>
00029 #include <OpenFOAM/addToRunTimeSelectionTable.H>
00030 #include <meshTools/triangleFuncs.H>
00031 #include <OpenFOAM/matchPoints.H>
00032 #include <OpenFOAM/globalIndex.H>
00033 #include <OpenFOAM/Time.H>
00034 
00035 #include <OpenFOAM/IFstream.H>
00036 #include <decompositionMethods/decompositionMethod.H>
00037 #include <OpenFOAM/vectorList.H>
00038 #include <OpenFOAM/PackedBoolList.H>
00039 
00040 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
00041 
00042 namespace Foam
00043 {
00044     defineTypeNameAndDebug(distributedTriSurfaceMesh, 0);
00045     addToRunTimeSelectionTable
00046     (
00047         searchableSurface,
00048         distributedTriSurfaceMesh,
00049         dict
00050     );
00051 }
00052 
00053 
00054 template<>
00055 const char*
00056 Foam::NamedEnum<Foam::distributedTriSurfaceMesh::distributionType, 3>::names[] =
00057 {
00058     "follow",
00059     "independent",
00060     "frozen"
00061 };
00062 
00063 const Foam::NamedEnum<Foam::distributedTriSurfaceMesh::distributionType, 3>
00064     Foam::distributedTriSurfaceMesh::distributionTypeNames_;
00065 
00066 
00067 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
00068 
00069 // Read my additional data from the dictionary
00070 bool Foam::distributedTriSurfaceMesh::read()
00071 {
00072     // Get bb of all domains.
00073     procBb_.setSize(Pstream::nProcs());
00074 
00075     procBb_[Pstream::myProcNo()] = List<treeBoundBox>(dict_.lookup("bounds"));
00076     Pstream::gatherList(procBb_);
00077     Pstream::scatterList(procBb_);
00078 
00079     // Distribution type
00080     distType_ = distributionTypeNames_.read(dict_.lookup("distributionType"));
00081 
00082     // Merge distance
00083     mergeDist_ = readScalar(dict_.lookup("mergeDistance"));
00084 
00085     return true;
00086 }
00087 
00088 
00089 // Is segment fully local?
00090 bool Foam::distributedTriSurfaceMesh::isLocal
00091 (
00092     const List<treeBoundBox>& myBbs,
00093     const point& start,
00094     const point& end
00095 )
00096 {
00097     forAll(myBbs, bbI)
00098     {
00099         if (myBbs[bbI].contains(start) && myBbs[bbI].contains(end))
00100         {
00101             return true;
00102         }
00103     }
00104     return false;
00105 }
00106 
00107 
00108 //void Foam::distributedTriSurfaceMesh::splitSegment
00109 //(
00110 //    const label segmentI,
00111 //    const point& start,
00112 //    const point& end,
00113 //    const treeBoundBox& bb,
00114 //
00115 //    DynamicList<segment>& allSegments,
00116 //    DynamicList<label>& allSegmentMap,
00117 //    DynamicList<label> sendMap
00118 //) const
00119 //{
00120 //    // Work points
00121 //    point clipPt0, clipPt1;
00122 //
00123 //    if (bb.contains(start))
00124 //    {
00125 //        // start within, trim end to bb
00126 //        bool clipped = bb.intersects(end, start, clipPt0);
00127 //
00128 //        if (clipped)
00129 //        {
00130 //            // segment from start to clippedStart passes
00131 //            // through proc.
00132 //            sendMap[procI].append(allSegments.size());
00133 //            allSegmentMap.append(segmentI);
00134 //            allSegments.append(segment(start, clipPt0));
00135 //        }
00136 //    }
00137 //    else if (bb.contains(end))
00138 //    {
00139 //        // end within, trim start to bb
00140 //        bool clipped = bb.intersects(start, end, clipPt0);
00141 //
00142 //        if (clipped)
00143 //        {
00144 //            sendMap[procI].append(allSegments.size());
00145 //            allSegmentMap.append(segmentI);
00146 //            allSegments.append(segment(clipPt0, end));
00147 //        }
00148 //    }
00149 //    else
00150 //    {
00151 //        // trim both
00152 //        bool clippedStart = bb.intersects(start, end, clipPt0);
00153 //
00154 //        if (clippedStart)
00155 //        {
00156 //            bool clippedEnd = bb.intersects(end, clipPt0, clipPt1);
00157 //
00158 //            if (clippedEnd)
00159 //            {
00160 //                // middle part of segment passes through proc.
00161 //                sendMap[procI].append(allSegments.size());
00162 //                allSegmentMap.append(segmentI);
00163 //                allSegments.append(segment(clipPt0, clipPt1));
00164 //            }
00165 //        }
00166 //    }
00167 //}
00168 
00169 
00170 void Foam::distributedTriSurfaceMesh::distributeSegment
00171 (
00172     const label segmentI,
00173     const point& start,
00174     const point& end,
00175 
00176     DynamicList<segment>& allSegments,
00177     DynamicList<label>& allSegmentMap,
00178     List<DynamicList<label> >& sendMap
00179 ) const
00180 {
00181     // Work points
00182     point clipPt;
00183 
00184 
00185     // 1. Fully local already handled outside. Note: retest is cheap.
00186     if (isLocal(procBb_[Pstream::myProcNo()], start, end))
00187     {
00188         return;
00189     }
00190 
00191 
00192     // 2. If fully inside one other processor, then only need to send
00193     // to that one processor even if it intersects another. Rare occurrence
00194     // but cheap to test.
00195     forAll(procBb_, procI)
00196     {
00197         if (procI != Pstream::myProcNo())
00198         {
00199             const List<treeBoundBox>& bbs = procBb_[procI];
00200 
00201             if (isLocal(bbs, start, end))
00202             {
00203                 sendMap[procI].append(allSegments.size());
00204                 allSegmentMap.append(segmentI);
00205                 allSegments.append(segment(start, end));
00206                 return;
00207             }
00208         }
00209     }
00210 
00211 
00212     // 3. If not contained in single processor send to all intersecting
00213     // processors.
00214     forAll(procBb_, procI)
00215     {
00216         const List<treeBoundBox>& bbs = procBb_[procI];
00217 
00218         forAll(bbs, bbI)
00219         {
00220             const treeBoundBox& bb = bbs[bbI];
00221 
00222             // Scheme a: any processor that intersects the segment gets
00223             // the segment.
00224 
00225             if (bb.intersects(start, end, clipPt))
00226             {
00227                 sendMap[procI].append(allSegments.size());
00228                 allSegmentMap.append(segmentI);
00229                 allSegments.append(segment(start, end));
00230             }
00231 
00232             // Alternative: any processor only gets clipped bit of
00233             // segment. This gives small problems with additional
00234             // truncation errors.
00235             //splitSegment
00236             //(
00237             //    segmentI,
00238             //    start,
00239             //    end,
00240             //    bb,
00241             //
00242             //    allSegments,
00243             //    allSegmentMap,
00244             //   sendMap[procI]
00245             //);
00246         }
00247     }
00248 }
00249 
00250 
00251 Foam::autoPtr<Foam::mapDistribute>
00252 Foam::distributedTriSurfaceMesh::distributeSegments
00253 (
00254     const pointField& start,
00255     const pointField& end,
00256 
00257     List<segment>& allSegments,
00258     labelList& allSegmentMap
00259 ) const
00260 {
00261     // Determine send map
00262     // ~~~~~~~~~~~~~~~~~~
00263 
00264     labelListList sendMap(Pstream::nProcs());
00265 
00266     {
00267         // Since intersection test is quite expensive compared to memory
00268         // allocation we use DynamicList to immediately store the segment
00269         // in the correct bin.
00270 
00271         // Segments to test
00272         DynamicList<segment> dynAllSegments(start.size());
00273         // Original index of segment
00274         DynamicList<label> dynAllSegmentMap(start.size());
00275         // Per processor indices into allSegments to send
00276         List<DynamicList<label> > dynSendMap(Pstream::nProcs());
00277 
00278         forAll(start, segmentI)
00279         {
00280             distributeSegment
00281             (
00282                 segmentI,
00283                 start[segmentI],
00284                 end[segmentI],
00285 
00286                 dynAllSegments,
00287                 dynAllSegmentMap,
00288                 dynSendMap
00289             );
00290         }
00291 
00292         // Convert dynamicList to labelList
00293         sendMap.setSize(Pstream::nProcs());
00294         forAll(sendMap, procI)
00295         {
00296             dynSendMap[procI].shrink();
00297             sendMap[procI].transfer(dynSendMap[procI]);
00298         }
00299 
00300         allSegments.transfer(dynAllSegments.shrink());
00301         allSegmentMap.transfer(dynAllSegmentMap.shrink());
00302     }
00303 
00304 
00305     // Send over how many I need to receive.
00306     labelListList sendSizes(Pstream::nProcs());
00307     sendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
00308     forAll(sendMap, procI)
00309     {
00310         sendSizes[Pstream::myProcNo()][procI] = sendMap[procI].size();
00311     }
00312     Pstream::gatherList(sendSizes);
00313     Pstream::scatterList(sendSizes);
00314 
00315 
00316     // Determine order of receiving
00317     labelListList constructMap(Pstream::nProcs());
00318 
00319     // My local segments first
00320     constructMap[Pstream::myProcNo()] = identity
00321     (
00322         sendMap[Pstream::myProcNo()].size()
00323     );
00324 
00325     label segmentI = constructMap[Pstream::myProcNo()].size();
00326     forAll(constructMap, procI)
00327     {
00328         if (procI != Pstream::myProcNo())
00329         {
00330             // What I need to receive is what other processor is sending to me.
00331             label nRecv = sendSizes[procI][Pstream::myProcNo()];
00332             constructMap[procI].setSize(nRecv);
00333 
00334             for (label i = 0; i < nRecv; i++)
00335             {
00336                 constructMap[procI][i] = segmentI++;
00337             }
00338         }
00339     }
00340 
00341     return autoPtr<mapDistribute>
00342     (
00343         new mapDistribute
00344         (
00345             segmentI,       // size after construction
00346             sendMap,
00347             constructMap,
00348             true            // reuse storage
00349         )
00350     );
00351 }
00352 
00353 
00354 void Foam::distributedTriSurfaceMesh::findLine
00355 (
00356     const bool nearestIntersection,
00357     const pointField& start,
00358     const pointField& end,
00359     List<pointIndexHit>& info
00360 ) const
00361 {
00362     const indexedOctree<treeDataTriSurface>& octree = tree();
00363 
00364     // Important:force synchronised construction of indexing
00365     const globalIndex& triIndexer = globalTris();
00366 
00367     // Initialise
00368     info.setSize(start.size());
00369     forAll(info, i)
00370     {
00371         info[i].setMiss();
00372     }
00373 
00374 
00375     // Do any local queries
00376     // ~~~~~~~~~~~~~~~~~~~~
00377 
00378     label nLocal = 0;
00379 
00380     forAll(start, i)
00381     {
00382         if (isLocal(procBb_[Pstream::myProcNo()], start[i], end[i]))
00383         {
00384             if (nearestIntersection)
00385             {
00386                 info[i] = octree.findLine(start[i], end[i]);
00387             }
00388             else
00389             {
00390                 info[i] = octree.findLineAny(start[i], end[i]);
00391             }
00392 
00393             if (info[i].hit())
00394             {
00395                 info[i].setIndex(triIndexer.toGlobal(info[i].index()));
00396             }
00397             nLocal++;
00398         }
00399     }
00400 
00401 
00402     if
00403     (
00404         Pstream::parRun()
00405      && (
00406             returnReduce(nLocal, sumOp<label>())
00407           < returnReduce(start.size(), sumOp<label>())
00408         )
00409     )
00410     {
00411         // Not all can be resolved locally. Build segments and map, send over
00412         // segments, do intersections, send back and merge.
00413 
00414 
00415         // Construct queries (segments)
00416         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00417 
00418         // Segments to test
00419         List<segment> allSegments(start.size());
00420         // Original index of segment
00421         labelList allSegmentMap(start.size());
00422 
00423         const autoPtr<mapDistribute> mapPtr
00424         (
00425             distributeSegments
00426             (
00427                 start,
00428                 end,
00429                 allSegments,
00430                 allSegmentMap
00431             )
00432         );
00433         const mapDistribute& map = mapPtr();
00434 
00435         label nOldAllSegments = allSegments.size();
00436 
00437 
00438         // Exchange the segments
00439         // ~~~~~~~~~~~~~~~~~~~~~
00440 
00441         map.distribute
00442         (
00443             Pstream::nonBlocking,   //Pstream::scheduled,
00444             List<labelPair>(0),     //map.schedule(),
00445             map.constructSize(),
00446             map.subMap(),           // what to send
00447             map.constructMap(),     // what to receive
00448             allSegments
00449         );
00450 
00451 
00452         // Do tests I need to do
00453         // ~~~~~~~~~~~~~~~~~~~~~
00454 
00455         // Intersections
00456         List<pointIndexHit> intersections(allSegments.size());
00457 
00458         forAll(allSegments, i)
00459         {
00460             if (nearestIntersection)
00461             {
00462                 intersections[i] = octree.findLine
00463                 (
00464                     allSegments[i].first(),
00465                     allSegments[i].second()
00466                 );
00467             }
00468             else
00469             {
00470                 intersections[i] = octree.findLineAny
00471                 (
00472                     allSegments[i].first(),
00473                     allSegments[i].second()
00474                 );
00475             }
00476 
00477             // Convert triangle index to global numbering
00478             if (intersections[i].hit())
00479             {
00480                 intersections[i].setIndex
00481                 (
00482                     triIndexer.toGlobal(intersections[i].index())
00483                 );
00484             }
00485         }
00486 
00487 
00488         // Exchange the intersections (opposite to segments)
00489         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00490 
00491         map.distribute
00492         (
00493             //Pstream::scheduled,
00494             //map.schedule            // Note reverse schedule
00495             //(
00496             //    map.constructMap(),
00497             //    map.subMap()
00498             //),
00499             Pstream::nonBlocking,
00500             List<labelPair>(0),
00501             nOldAllSegments,
00502             map.constructMap(),     // what to send
00503             map.subMap(),           // what to receive
00504             intersections
00505         );
00506 
00507 
00508         // Extract the hits
00509         // ~~~~~~~~~~~~~~~~
00510 
00511         forAll(intersections, i)
00512         {
00513             const pointIndexHit& allInfo = intersections[i];
00514             label segmentI = allSegmentMap[i];
00515             pointIndexHit& hitInfo = info[segmentI];
00516 
00517             if (allInfo.hit())
00518             {
00519                 if (!hitInfo.hit())
00520                 {
00521                     // No intersection yet so take this one
00522                     hitInfo = allInfo;
00523                 }
00524                 else if (nearestIntersection)
00525                 {
00526                     // Nearest intersection
00527                     if
00528                     (
00529                         magSqr(allInfo.hitPoint()-start[segmentI])
00530                       < magSqr(hitInfo.hitPoint()-start[segmentI])
00531                     )
00532                     {
00533                         hitInfo = allInfo;
00534                     }
00535                 }
00536             }
00537         }
00538     }
00539 }
00540 
00541 
00542 // Exchanges indices to the processor they come from.
00543 // - calculates exchange map
00544 // - uses map to calculate local triangle index
00545 Foam::autoPtr<Foam::mapDistribute>
00546 Foam::distributedTriSurfaceMesh::calcLocalQueries
00547 (
00548     const List<pointIndexHit>& info,
00549     labelList& triangleIndex
00550 ) const
00551 {
00552     triangleIndex.setSize(info.size());
00553 
00554     const globalIndex& triIndexer = globalTris();
00555 
00556 
00557     // Determine send map
00558     // ~~~~~~~~~~~~~~~~~~
00559 
00560     // Since determining which processor the query should go to is
00561     // cheap we do a multi-pass algorithm to save some memory temporarily.
00562 
00563     // 1. Count
00564     labelList nSend(Pstream::nProcs(), 0);
00565 
00566     forAll(info, i)
00567     {
00568         if (info[i].hit())
00569         {
00570             label procI = triIndexer.whichProcID(info[i].index());
00571             nSend[procI]++;
00572         }
00573     }
00574 
00575     // 2. Size sendMap
00576     labelListList sendMap(Pstream::nProcs());
00577     forAll(nSend, procI)
00578     {
00579         sendMap[procI].setSize(nSend[procI]);
00580         nSend[procI] = 0;
00581     }
00582 
00583     // 3. Fill sendMap
00584     forAll(info, i)
00585     {
00586         if (info[i].hit())
00587         {
00588             label procI = triIndexer.whichProcID(info[i].index());
00589             triangleIndex[i] = triIndexer.toLocal(procI, info[i].index());
00590             sendMap[procI][nSend[procI]++] = i;
00591         }
00592         else
00593         {
00594             triangleIndex[i] = -1;
00595         }
00596     }
00597 
00598 
00599     // Send over how many I need to receive
00600     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00601 
00602     labelListList sendSizes(Pstream::nProcs());
00603     sendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
00604     forAll(sendMap, procI)
00605     {
00606         sendSizes[Pstream::myProcNo()][procI] = sendMap[procI].size();
00607     }
00608     Pstream::gatherList(sendSizes);
00609     Pstream::scatterList(sendSizes);
00610 
00611 
00612     // Determine receive map
00613     // ~~~~~~~~~~~~~~~~~~~~~
00614 
00615     labelListList constructMap(Pstream::nProcs());
00616 
00617     // My local segments first
00618     constructMap[Pstream::myProcNo()] = identity
00619     (
00620         sendMap[Pstream::myProcNo()].size()
00621     );
00622 
00623     label segmentI = constructMap[Pstream::myProcNo()].size();
00624     forAll(constructMap, procI)
00625     {
00626         if (procI != Pstream::myProcNo())
00627         {
00628             // What I need to receive is what other processor is sending to me.
00629             label nRecv = sendSizes[procI][Pstream::myProcNo()];
00630             constructMap[procI].setSize(nRecv);
00631 
00632             for (label i = 0; i < nRecv; i++)
00633             {
00634                 constructMap[procI][i] = segmentI++;
00635             }
00636         }
00637     }
00638 
00639 
00640     // Pack into distribution map
00641     // ~~~~~~~~~~~~~~~~~~~~~~~~~~
00642 
00643     autoPtr<mapDistribute> mapPtr
00644     (
00645         new mapDistribute
00646         (
00647             segmentI,       // size after construction
00648             sendMap,
00649             constructMap,
00650             true            // reuse storage
00651         )
00652     );
00653     const mapDistribute& map = mapPtr();
00654 
00655 
00656     // Send over queries
00657     // ~~~~~~~~~~~~~~~~~
00658 
00659     map.distribute
00660     (
00661         //Pstream::scheduled,
00662         //map.schedule(),
00663         Pstream::nonBlocking,
00664         List<labelPair>(0),
00665         map.constructSize(),
00666         map.subMap(),           // what to send
00667         map.constructMap(),     // what to receive
00668         triangleIndex
00669     );
00670 
00671 
00672     return mapPtr;
00673 }
00674 
00675 
00676 Foam::label Foam::distributedTriSurfaceMesh::calcOverlappingProcs
00677 (
00678     const point& centre,
00679     const scalar radiusSqr,
00680     boolList& overlaps
00681 ) const
00682 {
00683     overlaps = false;
00684     label nOverlaps = 0;
00685 
00686     forAll(procBb_, procI)
00687     {
00688         const List<treeBoundBox>& bbs = procBb_[procI];
00689 
00690         forAll(bbs, bbI)
00691         {
00692             if (bbs[bbI].overlaps(centre, radiusSqr))
00693             {
00694                 overlaps[procI] = true;
00695                 nOverlaps++;
00696                 break;
00697             }
00698         }
00699     }
00700     return nOverlaps;
00701 }
00702 
00703 
00704 // Generate queries for parallel distance calculation
00705 // - calculates exchange map
00706 // - uses map to exchange points and radius
00707 Foam::autoPtr<Foam::mapDistribute>
00708 Foam::distributedTriSurfaceMesh::calcLocalQueries
00709 (
00710     const pointField& centres,
00711     const scalarField& radiusSqr,
00712 
00713     pointField& allCentres,
00714     scalarField& allRadiusSqr,
00715     labelList& allSegmentMap
00716 ) const
00717 {
00718     // Determine queries
00719     // ~~~~~~~~~~~~~~~~~
00720 
00721     labelListList sendMap(Pstream::nProcs());
00722 
00723     {
00724         // Queries
00725         DynamicList<point> dynAllCentres(centres.size());
00726         DynamicList<scalar> dynAllRadiusSqr(centres.size());
00727         // Original index of segment
00728         DynamicList<label> dynAllSegmentMap(centres.size());
00729         // Per processor indices into allSegments to send
00730         List<DynamicList<label> > dynSendMap(Pstream::nProcs());
00731 
00732         // Work array - whether processor bb overlaps the bounding sphere.
00733         boolList procBbOverlaps(Pstream::nProcs());
00734 
00735         forAll(centres, centreI)
00736         {
00737             // Find the processor this sample+radius overlaps.
00738             calcOverlappingProcs
00739             (
00740                 centres[centreI],
00741                 radiusSqr[centreI],
00742                 procBbOverlaps
00743             );
00744 
00745             forAll(procBbOverlaps, procI)
00746             {
00747                 if (procI != Pstream::myProcNo() && procBbOverlaps[procI])
00748                 {
00749                     dynSendMap[procI].append(dynAllCentres.size());
00750                     dynAllSegmentMap.append(centreI);
00751                     dynAllCentres.append(centres[centreI]);
00752                     dynAllRadiusSqr.append(radiusSqr[centreI]);
00753                 }
00754             }
00755         }
00756 
00757         // Convert dynamicList to labelList
00758         sendMap.setSize(Pstream::nProcs());
00759         forAll(sendMap, procI)
00760         {
00761             dynSendMap[procI].shrink();
00762             sendMap[procI].transfer(dynSendMap[procI]);
00763         }
00764 
00765         allCentres.transfer(dynAllCentres.shrink());
00766         allRadiusSqr.transfer(dynAllRadiusSqr.shrink());
00767         allSegmentMap.transfer(dynAllSegmentMap.shrink());
00768     }
00769 
00770 
00771     // Send over how many I need to receive.
00772     labelListList sendSizes(Pstream::nProcs());
00773     sendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
00774     forAll(sendMap, procI)
00775     {
00776         sendSizes[Pstream::myProcNo()][procI] = sendMap[procI].size();
00777     }
00778     Pstream::gatherList(sendSizes);
00779     Pstream::scatterList(sendSizes);
00780 
00781 
00782     // Determine order of receiving
00783     labelListList constructMap(Pstream::nProcs());
00784 
00785     // My local segments first
00786     constructMap[Pstream::myProcNo()] = identity
00787     (
00788         sendMap[Pstream::myProcNo()].size()
00789     );
00790 
00791     label segmentI = constructMap[Pstream::myProcNo()].size();
00792     forAll(constructMap, procI)
00793     {
00794         if (procI != Pstream::myProcNo())
00795         {
00796             // What I need to receive is what other processor is sending to me.
00797             label nRecv = sendSizes[procI][Pstream::myProcNo()];
00798             constructMap[procI].setSize(nRecv);
00799 
00800             for (label i = 0; i < nRecv; i++)
00801             {
00802                 constructMap[procI][i] = segmentI++;
00803             }
00804         }
00805     }
00806 
00807     autoPtr<mapDistribute> mapPtr
00808     (
00809         new mapDistribute
00810         (
00811             segmentI,       // size after construction
00812             sendMap,
00813             constructMap,
00814             true            // reuse storage
00815         )
00816     );
00817     return mapPtr;
00818 }
00819 
00820 
00821 // Find bounding boxes that guarantee a more or less uniform distribution
00822 // of triangles. Decomposition in here is only used to get the bounding
00823 // boxes, actual decomposition is done later on.
00824 // Returns a per processor a list of bounding boxes that most accurately
00825 // describe the shape. For now just a single bounding box per processor but
00826 // optimisation might be to determine a better fitting shape.
00827 Foam::List<Foam::List<Foam::treeBoundBox> >
00828 Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
00829 (
00830     const triSurface& s
00831 )
00832 {
00833     if (!decomposer_.valid())
00834     {
00835         // Use current decomposer.
00836         // Note: or always use hierarchical?
00837         IOdictionary decomposeDict
00838         (
00839             IOobject
00840             (
00841                 "decomposeParDict",
00842                 searchableSurface::time().system(),
00843                 searchableSurface::time(),
00844                 IOobject::MUST_READ,
00845                 IOobject::NO_WRITE,
00846                 false
00847             )
00848         );
00849         decomposer_ = decompositionMethod::New(decomposeDict);
00850 
00851         if (!decomposer_().parallelAware())
00852         {
00853             FatalErrorIn
00854             (
00855                 "distributedTriSurfaceMesh::independentlyDistributedBbs"
00856                 "(const triSurface&)"
00857             )   << "The decomposition method " << decomposer_().typeName
00858                 << " does not decompose in parallel."
00859                 << " Please choose one that does." << exit(FatalError);
00860         }
00861     }
00862 
00863     // Do decomposition according to triangle centre
00864     pointField triCentres(s.size());
00865     forAll (s, triI)
00866     {
00867         triCentres[triI] = s[triI].centre(s.points());
00868     }
00869 
00870     // Do the actual decomposition
00871     labelList distribution(decomposer_->decompose(triCentres));
00872 
00873     // Find bounding box for all triangles on new distribution.
00874 
00875     // Initialise to inverted box (VGREAT, -VGREAT)
00876     List<List<treeBoundBox> > bbs(Pstream::nProcs());
00877     forAll(bbs, procI)
00878     {
00879         bbs[procI].setSize(1);
00880         //bbs[procI][0] = boundBox::invertedBox;
00881         bbs[procI][0].min() = point( VGREAT,  VGREAT,  VGREAT);
00882         bbs[procI][0].max() = point(-VGREAT, -VGREAT, -VGREAT);
00883     }
00884 
00885     forAll (s, triI)
00886     {
00887         point& bbMin = bbs[distribution[triI]][0].min();
00888         point& bbMax = bbs[distribution[triI]][0].max();
00889 
00890         const labelledTri& f = s[triI];
00891         const point& p0 = s.points()[f[0]];
00892         const point& p1 = s.points()[f[1]];
00893         const point& p2 = s.points()[f[2]];
00894 
00895         bbMin = min(bbMin, p0);
00896         bbMin = min(bbMin, p1);
00897         bbMin = min(bbMin, p2);
00898 
00899         bbMax = max(bbMax, p0);
00900         bbMax = max(bbMax, p1);
00901         bbMax = max(bbMax, p2);
00902     }
00903 
00904     // Now combine for all processors and convert to correct format.
00905     forAll(bbs, procI)
00906     {
00907         forAll(bbs[procI], i)
00908         {
00909             reduce(bbs[procI][i].min(), minOp<point>());
00910             reduce(bbs[procI][i].max(), maxOp<point>());
00911         }
00912     }
00913     return bbs;
00914 }
00915 
00916 
00917 // Does any part of triangle overlap bb.
00918 bool Foam::distributedTriSurfaceMesh::overlaps
00919 (
00920     const List<treeBoundBox>& bbs,
00921     const point& p0,
00922     const point& p1,
00923     const point& p2
00924 )
00925 {
00926     forAll(bbs, bbI)
00927     {
00928         const treeBoundBox& bb = bbs[bbI];
00929 
00930         treeBoundBox triBb(p0, p0);
00931         triBb.min() = min(triBb.min(), p1);
00932         triBb.min() = min(triBb.min(), p2);
00933 
00934         triBb.max() = max(triBb.max(), p1);
00935         triBb.max() = max(triBb.max(), p2);
00936 
00937         //- Exact test of triangle intersecting bb
00938 
00939         // Quick rejection. If whole bounding box of tri is outside cubeBb then
00940         // there will be no intersection.
00941         if (bb.overlaps(triBb))
00942         {
00943             // Check if one or more triangle point inside
00944             if (bb.contains(p0) || bb.contains(p1) || bb.contains(p2))
00945             {
00946                 // One or more points inside
00947                 return true;
00948             }
00949 
00950             // Now we have the difficult case: all points are outside but
00951             // connecting edges might go through cube. Use fast intersection
00952             // of bounding box.
00953             bool intersect = triangleFuncs::intersectBb(p0, p1, p2, bb);
00954 
00955             if (intersect)
00956             {
00957                 return true;
00958             }
00959         }
00960     }
00961     return false;
00962 }
00963 
00964 
00965 void Foam::distributedTriSurfaceMesh::subsetMeshMap
00966 (
00967     const triSurface& s,
00968     const boolList& include,
00969     const label nIncluded,
00970     labelList& newToOldPoints,
00971     labelList& oldToNewPoints,
00972     labelList& newToOldFaces
00973 )
00974 {
00975     newToOldFaces.setSize(nIncluded);
00976     newToOldPoints.setSize(s.points().size());
00977     oldToNewPoints.setSize(s.points().size());
00978     oldToNewPoints = -1;
00979     {
00980         label faceI = 0;
00981         label pointI = 0;
00982 
00983         forAll(include, oldFacei)
00984         {
00985             if (include[oldFacei])
00986             {
00987                 // Store new faces compact
00988                 newToOldFaces[faceI++] = oldFacei;
00989 
00990                 // Renumber labels for triangle
00991                 const labelledTri& tri = s[oldFacei];
00992 
00993                 forAll(tri, fp)
00994                 {
00995                     label oldPointI = tri[fp];
00996 
00997                     if (oldToNewPoints[oldPointI] == -1)
00998                     {
00999                         oldToNewPoints[oldPointI] = pointI;
01000                         newToOldPoints[pointI++] = oldPointI;
01001                     }
01002                 }
01003             }
01004         }
01005         newToOldPoints.setSize(pointI);
01006     }
01007 }
01008 
01009 
01010 Foam::triSurface Foam::distributedTriSurfaceMesh::subsetMesh
01011 (
01012     const triSurface& s,
01013     const labelList& newToOldPoints,
01014     const labelList& oldToNewPoints,
01015     const labelList& newToOldFaces
01016 )
01017 {
01018     // Extract points
01019     pointField newPoints(newToOldPoints.size());
01020     forAll(newToOldPoints, i)
01021     {
01022         newPoints[i] = s.points()[newToOldPoints[i]];
01023     }
01024     // Extract faces
01025     List<labelledTri> newTriangles(newToOldFaces.size());
01026 
01027     forAll(newToOldFaces, i)
01028     {
01029         // Get old vertex labels
01030         const labelledTri& tri = s[newToOldFaces[i]];
01031 
01032         newTriangles[i][0] = oldToNewPoints[tri[0]];
01033         newTriangles[i][1] = oldToNewPoints[tri[1]];
01034         newTriangles[i][2] = oldToNewPoints[tri[2]];
01035         newTriangles[i].region() = tri.region();
01036     }
01037 
01038     // Reuse storage.
01039     return triSurface(newTriangles, s.patches(), newPoints, true);
01040 }
01041 
01042 
01043 Foam::triSurface Foam::distributedTriSurfaceMesh::subsetMesh
01044 (
01045     const triSurface& s,
01046     const boolList& include,
01047     labelList& newToOldPoints,
01048     labelList& newToOldFaces
01049 )
01050 {
01051     label n = 0;
01052 
01053     forAll(include, i)
01054     {
01055         if (include[i])
01056         {
01057             n++;
01058         }
01059     }
01060 
01061     labelList oldToNewPoints;
01062     subsetMeshMap
01063     (
01064         s,
01065         include,
01066         n,
01067         newToOldPoints,
01068         oldToNewPoints,
01069         newToOldFaces
01070     );
01071 
01072     return subsetMesh
01073     (
01074         s,
01075         newToOldPoints,
01076         oldToNewPoints,
01077         newToOldFaces
01078     );
01079 }
01080 
01081 
01082 Foam::triSurface Foam::distributedTriSurfaceMesh::subsetMesh
01083 (
01084     const triSurface& s,
01085     const labelList& newToOldFaces,
01086     labelList& newToOldPoints
01087 )
01088 {
01089     const boolList include
01090     (
01091         createWithValues<boolList>
01092         (
01093             s.size(),
01094             false,
01095             newToOldFaces,
01096             true
01097         )
01098     );
01099 
01100     newToOldPoints.setSize(s.points().size());
01101     labelList oldToNewPoints(s.points().size(), -1);
01102     {
01103         label pointI = 0;
01104 
01105         forAll(include, oldFacei)
01106         {
01107             if (include[oldFacei])
01108             {
01109                 // Renumber labels for triangle
01110                 const labelledTri& tri = s[oldFacei];
01111 
01112                 forAll(tri, fp)
01113                 {
01114                     label oldPointI = tri[fp];
01115 
01116                     if (oldToNewPoints[oldPointI] == -1)
01117                     {
01118                         oldToNewPoints[oldPointI] = pointI;
01119                         newToOldPoints[pointI++] = oldPointI;
01120                     }
01121                 }
01122             }
01123         }
01124         newToOldPoints.setSize(pointI);
01125     }
01126 
01127     return subsetMesh
01128     (
01129         s,
01130         newToOldPoints,
01131         oldToNewPoints,
01132         newToOldFaces
01133     );
01134 }
01135 
01136 
01137 Foam::label Foam::distributedTriSurfaceMesh::findTriangle
01138 (
01139     const List<labelledTri>& allFaces,
01140     const labelListList& allPointFaces,
01141     const labelledTri& otherF
01142 )
01143 {
01144     // allFaces connected to otherF[0]
01145     const labelList& pFaces = allPointFaces[otherF[0]];
01146 
01147     forAll(pFaces, i)
01148     {
01149         const labelledTri& f = allFaces[pFaces[i]];
01150 
01151         if (f.region() == otherF.region())
01152         {
01153             // Find index of otherF[0]
01154             label fp0 = findIndex(f, otherF[0]);
01155             // Check rest of triangle in same order
01156             label fp1 = f.fcIndex(fp0);
01157             label fp2 = f.fcIndex(fp1);
01158 
01159             if (f[fp1] == otherF[1] && f[fp2] == otherF[2])
01160             {
01161                 return pFaces[i];
01162             }
01163         }
01164     }
01165     return -1;
01166 }
01167 
01168 
01169 // Merge into allSurf.
01170 void Foam::distributedTriSurfaceMesh::merge
01171 (
01172     const scalar mergeDist,
01173     const List<labelledTri>& subTris,
01174     const pointField& subPoints,
01175 
01176     List<labelledTri>& allTris,
01177     pointField& allPoints,
01178 
01179     labelList& faceConstructMap,
01180     labelList& pointConstructMap
01181 )
01182 {
01183     labelList subToAll;
01184     matchPoints
01185     (
01186         subPoints,
01187         allPoints,
01188         scalarField(subPoints.size(), mergeDist),   // match distance
01189         false,                                      // verbose
01190         pointConstructMap
01191     );
01192 
01193     label nOldAllPoints = allPoints.size();
01194 
01195 
01196     // Add all unmatched points
01197     // ~~~~~~~~~~~~~~~~~~~~~~~~
01198 
01199     label allPointI = nOldAllPoints;
01200     forAll(pointConstructMap, pointI)
01201     {
01202         if (pointConstructMap[pointI] == -1)
01203         {
01204             pointConstructMap[pointI] = allPointI++;
01205         }
01206     }
01207 
01208     if (allPointI > nOldAllPoints)
01209     {
01210         allPoints.setSize(allPointI);
01211 
01212         forAll(pointConstructMap, pointI)
01213         {
01214             if (pointConstructMap[pointI] >= nOldAllPoints)
01215             {
01216                 allPoints[pointConstructMap[pointI]] = subPoints[pointI];
01217             }
01218         }
01219     }
01220 
01221 
01222     // To check whether triangles are same we use pointFaces.
01223     labelListList allPointFaces;
01224     invertManyToMany(nOldAllPoints, allTris, allPointFaces);
01225 
01226 
01227     // Add all unmatched triangles
01228     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
01229 
01230     label allTriI = allTris.size();
01231     allTris.setSize(allTriI + subTris.size());
01232 
01233     faceConstructMap.setSize(subTris.size());
01234 
01235     forAll(subTris, triI)
01236     {
01237         const labelledTri& subTri = subTris[triI];
01238 
01239         // Get triangle in new numbering
01240         labelledTri mappedTri
01241         (
01242             pointConstructMap[subTri[0]],
01243             pointConstructMap[subTri[1]],
01244             pointConstructMap[subTri[2]],
01245             subTri.region()
01246         );
01247 
01248 
01249         // Check if all points were already in surface
01250         bool fullMatch = true;
01251 
01252         forAll(mappedTri, fp)
01253         {
01254             if (mappedTri[fp] >= nOldAllPoints)
01255             {
01256                 fullMatch = false;
01257                 break;
01258             }
01259         }
01260 
01261         if (fullMatch)
01262         {
01263             // All three points are mapped to old points. See if same
01264             // triangle.
01265             label i = findTriangle
01266             (
01267                 allTris,
01268                 allPointFaces,
01269                 mappedTri
01270             );
01271 
01272             if (i == -1)
01273             {
01274                 // Add
01275                 faceConstructMap[triI] = allTriI;
01276                 allTris[allTriI] = mappedTri;
01277                 allTriI++;
01278             }
01279             else
01280             {
01281                 faceConstructMap[triI] = i;
01282             }
01283         }
01284         else
01285         {
01286             // Add
01287             faceConstructMap[triI] = allTriI;
01288             allTris[allTriI] = mappedTri;
01289             allTriI++;
01290         }
01291     }
01292     allTris.setSize(allTriI);
01293 }
01294 
01295 
01296 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
01297 
01298 Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh
01299 (
01300     const IOobject& io,
01301     const triSurface& s,
01302     const dictionary& dict
01303 )
01304 :
01305     triSurfaceMesh(io, s),
01306     dict_
01307     (
01308         IOobject
01309         (
01310             searchableSurface::name() + "Dict",
01311             searchableSurface::instance(),
01312             searchableSurface::local(),
01313             searchableSurface::db(),
01314             searchableSurface::NO_READ,
01315             searchableSurface::writeOpt(),
01316             searchableSurface::registerObject()
01317         ),
01318         dict
01319     )
01320 {
01321     read();
01322 
01323     if (debug)
01324     {
01325         Info<< "Constructed from triSurface:" << endl;
01326         writeStats(Info);
01327 
01328         labelList nTris(Pstream::nProcs());
01329         nTris[Pstream::myProcNo()] = triSurface::size();
01330         Pstream::gatherList(nTris);
01331         Pstream::scatterList(nTris);
01332 
01333         Info<< endl<< "\tproc\ttris\tbb" << endl;
01334         forAll(nTris, procI)
01335         {
01336             Info<< '\t' << procI << '\t' << nTris[procI]
01337                  << '\t' << procBb_[procI] << endl;
01338         }
01339         Info<< endl;
01340     }
01341 }
01342 
01343 
01344 Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh(const IOobject& io)
01345 :
01346     //triSurfaceMesh(io),
01347     triSurfaceMesh
01348     (
01349         IOobject
01350         (
01351             io.name(),
01352             io.time().findInstance(io.local(), word::null),
01353             io.local(),
01354             io.db(),
01355             io.readOpt(),
01356             io.writeOpt(),
01357             io.registerObject()
01358         )
01359     ),
01360     dict_
01361     (
01362         IOobject
01363         (
01364             searchableSurface::name() + "Dict",
01365             searchableSurface::instance(),
01366             searchableSurface::local(),
01367             searchableSurface::db(),
01368             searchableSurface::readOpt(),
01369             searchableSurface::writeOpt(),
01370             searchableSurface::registerObject()
01371         )
01372     )
01373 {
01374     read();
01375 
01376     if (debug)
01377     {
01378         Info<< "Read distributedTriSurface from " << io.objectPath()
01379             << ':' << endl;
01380         writeStats(Info);
01381 
01382         labelList nTris(Pstream::nProcs());
01383         nTris[Pstream::myProcNo()] = triSurface::size();
01384         Pstream::gatherList(nTris);
01385         Pstream::scatterList(nTris);
01386 
01387         Info<< endl<< "\tproc\ttris\tbb" << endl;
01388         forAll(nTris, procI)
01389         {
01390             Info<< '\t' << procI << '\t' << nTris[procI]
01391                  << '\t' << procBb_[procI] << endl;
01392         }
01393         Info<< endl;
01394     }
01395 }
01396 
01397 
01398 Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh
01399 (
01400     const IOobject& io,
01401     const dictionary& dict
01402 )
01403 :
01404     //triSurfaceMesh(io, dict),
01405     triSurfaceMesh
01406     (
01407         IOobject
01408         (
01409             io.name(),
01410             io.time().findInstance(io.local(), word::null),
01411             io.local(),
01412             io.db(),
01413             io.readOpt(),
01414             io.writeOpt(),
01415             io.registerObject()
01416         ),
01417         dict
01418     ),
01419     dict_
01420     (
01421         IOobject
01422         (
01423             searchableSurface::name() + "Dict",
01424             searchableSurface::instance(),
01425             searchableSurface::local(),
01426             searchableSurface::db(),
01427             searchableSurface::readOpt(),
01428             searchableSurface::writeOpt(),
01429             searchableSurface::registerObject()
01430         )
01431     )
01432 {
01433     read();
01434 
01435     if (debug)
01436     {
01437         Info<< "Read distributedTriSurface from " << io.objectPath()
01438             << " and dictionary:" << endl;
01439         writeStats(Info);
01440 
01441         labelList nTris(Pstream::nProcs());
01442         nTris[Pstream::myProcNo()] = triSurface::size();
01443         Pstream::gatherList(nTris);
01444         Pstream::scatterList(nTris);
01445 
01446         Info<< endl<< "\tproc\ttris\tbb" << endl;
01447         forAll(nTris, procI)
01448         {
01449             Info<< '\t' << procI << '\t' << nTris[procI]
01450                  << '\t' << procBb_[procI] << endl;
01451         }
01452         Info<< endl;
01453     }
01454 }
01455 
01456 
01457 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
01458 
01459 Foam::distributedTriSurfaceMesh::~distributedTriSurfaceMesh()
01460 {
01461     clearOut();
01462 }
01463 
01464 
01465 void Foam::distributedTriSurfaceMesh::clearOut()
01466 {
01467     globalTris_.clear();
01468     triSurfaceMesh::clearOut();
01469 }
01470 
01471 
01472 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
01473 
01474 const Foam::globalIndex& Foam::distributedTriSurfaceMesh::globalTris() const
01475 {
01476     if (!globalTris_.valid())
01477     {
01478         globalTris_.reset(new globalIndex(triSurface::size()));
01479     }
01480     return globalTris_;
01481 }
01482 
01483 
01484 void Foam::distributedTriSurfaceMesh::findNearest
01485 (
01486     const pointField& samples,
01487     const scalarField& nearestDistSqr,
01488     List<pointIndexHit>& info
01489 ) const
01490 {
01491     const indexedOctree<treeDataTriSurface>& octree = tree();
01492 
01493     // Important:force synchronised construction of indexing
01494     const globalIndex& triIndexer = globalTris();
01495 
01496 
01497     // Initialise
01498     // ~~~~~~~~~~
01499 
01500     info.setSize(samples.size());
01501     forAll(info, i)
01502     {
01503         info[i].setMiss();
01504     }
01505 
01506 
01507 
01508     // Do any local queries
01509     // ~~~~~~~~~~~~~~~~~~~~
01510 
01511     label nLocal = 0;
01512 
01513     {
01514         // Work array - whether processor bb overlaps the bounding sphere.
01515         boolList procBbOverlaps(Pstream::nProcs());
01516 
01517         forAll(samples, i)
01518         {
01519             // Find the processor this sample+radius overlaps.
01520             label nProcs = calcOverlappingProcs
01521             (
01522                 samples[i],
01523                 nearestDistSqr[i],
01524                 procBbOverlaps
01525             );
01526 
01527             // Overlaps local processor?
01528             if (procBbOverlaps[Pstream::myProcNo()])
01529             {
01530                 info[i] = octree.findNearest(samples[i], nearestDistSqr[i]);
01531                 if (info[i].hit())
01532                 {
01533                     info[i].setIndex(triIndexer.toGlobal(info[i].index()));
01534                 }
01535                 if (nProcs == 1)
01536                 {
01537                     // Fully local
01538                     nLocal++;
01539                 }
01540             }
01541         }
01542     }
01543 
01544 
01545     if
01546     (
01547         Pstream::parRun()
01548      && (
01549             returnReduce(nLocal, sumOp<label>())
01550           < returnReduce(samples.size(), sumOp<label>())
01551         )
01552     )
01553     {
01554         // Not all can be resolved locally. Build queries and map, send over
01555         // queries, do intersections, send back and merge.
01556 
01557         // Calculate queries and exchange map
01558         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
01559 
01560         pointField allCentres;
01561         scalarField allRadiusSqr;
01562         labelList allSegmentMap;
01563         autoPtr<mapDistribute> mapPtr
01564         (
01565             calcLocalQueries
01566             (
01567                 samples,
01568                 nearestDistSqr,
01569 
01570                 allCentres,
01571                 allRadiusSqr,
01572                 allSegmentMap
01573             )
01574         );
01575         const mapDistribute& map = mapPtr();
01576 
01577 
01578         // swap samples to local processor
01579         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
01580 
01581         map.distribute
01582         (
01583             //Pstream::scheduled,
01584             //map.schedule(),
01585             Pstream::nonBlocking,
01586             List<labelPair>(0),
01587             map.constructSize(),
01588             map.subMap(),           // what to send
01589             map.constructMap(),     // what to receive
01590             allCentres
01591         );
01592         map.distribute
01593         (
01594             //Pstream::scheduled,
01595             //map.schedule(),
01596             Pstream::nonBlocking,
01597             List<labelPair>(0),
01598             map.constructSize(),
01599             map.subMap(),           // what to send
01600             map.constructMap(),     // what to receive
01601             allRadiusSqr
01602         );
01603 
01604 
01605         // Do my tests
01606         // ~~~~~~~~~~~
01607 
01608         List<pointIndexHit> allInfo(allCentres.size());
01609         forAll(allInfo, i)
01610         {
01611             allInfo[i] = octree.findNearest
01612             (
01613                 allCentres[i],
01614                 allRadiusSqr[i]
01615             );
01616             if (allInfo[i].hit())
01617             {
01618                 allInfo[i].setIndex(triIndexer.toGlobal(allInfo[i].index()));
01619             }
01620         }
01621 
01622 
01623         // Send back results
01624         // ~~~~~~~~~~~~~~~~~
01625 
01626         map.distribute
01627         (
01628             //Pstream::scheduled,
01629             //map.schedule            // note reverse schedule
01630             //(
01631             //    map.constructMap(),
01632             //    map.subMap()
01633             //),
01634             Pstream::nonBlocking,
01635             List<labelPair>(0),
01636             allSegmentMap.size(),
01637             map.constructMap(),     // what to send
01638             map.subMap(),           // what to receive
01639             allInfo
01640         );
01641 
01642 
01643         // Extract information
01644         // ~~~~~~~~~~~~~~~~~~~
01645 
01646         forAll(allInfo, i)
01647         {
01648             if (allInfo[i].hit())
01649             {
01650                 label pointI = allSegmentMap[i];
01651 
01652                 if (!info[pointI].hit())
01653                 {
01654                     // No intersection yet so take this one
01655                     info[pointI] = allInfo[i];
01656                 }
01657                 else
01658                 {
01659                     // Nearest intersection
01660                     if
01661                     (
01662                         magSqr(allInfo[i].hitPoint()-samples[pointI])
01663                       < magSqr(info[pointI].hitPoint()-samples[pointI])
01664                     )
01665                     {
01666                         info[pointI] = allInfo[i];
01667                     }
01668                 }
01669             }
01670         }
01671     }
01672 }
01673 
01674 
01675 void Foam::distributedTriSurfaceMesh::findLine
01676 (
01677     const pointField& start,
01678     const pointField& end,
01679     List<pointIndexHit>& info
01680 ) const
01681 {
01682     findLine
01683     (
01684         true,   // nearestIntersection
01685         start,
01686         end,
01687         info
01688     );
01689 }
01690 
01691 
01692 void Foam::distributedTriSurfaceMesh::findLineAny
01693 (
01694     const pointField& start,
01695     const pointField& end,
01696     List<pointIndexHit>& info
01697 ) const
01698 {
01699     findLine
01700     (
01701         true,   // nearestIntersection
01702         start,
01703         end,
01704         info
01705     );
01706 }
01707 
01708 
01709 void Foam::distributedTriSurfaceMesh::findLineAll
01710 (
01711     const pointField& start,
01712     const pointField& end,
01713     List<List<pointIndexHit> >& info
01714 ) const
01715 {
01716     // Reuse fineLine. We could modify all of findLine to do multiple
01717     // intersections but this would mean a lot of data transferred so
01718     // for now we just find nearest intersection and retest from that
01719     // intersection onwards.
01720 
01721     // Work array.
01722     List<pointIndexHit> hitInfo(start.size());
01723 
01724     findLine
01725     (
01726         true,   // nearestIntersection
01727         start,
01728         end,
01729         hitInfo
01730     );
01731 
01732     // Tolerances:
01733     // To find all intersections we add a small vector to the last intersection
01734     // This is chosen such that
01735     // - it is significant (SMALL is smallest representative relative tolerance;
01736     //   we need something bigger since we're doing calculations)
01737     // - if the start-end vector is zero we still progress
01738     const vectorField dirVec(end-start);
01739     const scalarField magSqrDirVec(magSqr(dirVec));
01740     const vectorField smallVec
01741     (
01742         Foam::sqrt(SMALL)*dirVec
01743       + vector(ROOTVSMALL,ROOTVSMALL,ROOTVSMALL)
01744     );
01745 
01746     // Copy to input and compact any hits
01747     labelList pointMap(start.size());
01748     pointField e0(start.size());
01749     pointField e1(start.size());
01750     label compactI = 0;
01751 
01752     info.setSize(hitInfo.size());
01753     forAll(hitInfo, pointI)
01754     {
01755         if (hitInfo[pointI].hit())
01756         {
01757             info[pointI].setSize(1);
01758             info[pointI][0] = hitInfo[pointI];
01759 
01760             point pt = hitInfo[pointI].hitPoint() + smallVec[pointI];
01761 
01762             if (((pt-start[pointI])&dirVec[pointI]) <= magSqrDirVec[pointI])
01763             {
01764                 e0[compactI] = pt;
01765                 e1[compactI] = end[pointI];
01766                 pointMap[compactI] = pointI;
01767                 compactI++;
01768             }
01769         }
01770         else
01771         {
01772             info[pointI].clear();
01773         }
01774     }
01775 
01776     e0.setSize(compactI);
01777     e1.setSize(compactI);
01778     pointMap.setSize(compactI);
01779 
01780     while (returnReduce(e0.size(), sumOp<label>()) > 0)
01781     {
01782         findLine
01783         (
01784             true,   // nearestIntersection
01785             e0,
01786             e1,
01787             hitInfo
01788         );
01789 
01790 
01791         // Extract
01792         label compactI = 0;
01793         forAll(hitInfo, i)
01794         {
01795             if (hitInfo[i].hit())
01796             {
01797                 label pointI = pointMap[i];
01798 
01799                 label sz = info[pointI].size();
01800                 info[pointI].setSize(sz+1);
01801                 info[pointI][sz] = hitInfo[i];
01802 
01803                 point pt = hitInfo[i].hitPoint() + smallVec[pointI];
01804 
01805                 if (((pt-start[pointI])&dirVec[pointI]) <= magSqrDirVec[pointI])
01806                 {
01807                     e0[compactI] = pt;
01808                     e1[compactI] = end[pointI];
01809                     pointMap[compactI] = pointI;
01810                     compactI++;
01811                 }
01812             }
01813         }
01814 
01815         // Trim
01816         e0.setSize(compactI);
01817         e1.setSize(compactI);
01818         pointMap.setSize(compactI);
01819     }
01820 }
01821 
01822 
01823 void Foam::distributedTriSurfaceMesh::getRegion
01824 (
01825     const List<pointIndexHit>& info,
01826     labelList& region
01827 ) const
01828 {
01829     if (!Pstream::parRun())
01830     {
01831         region.setSize(info.size());
01832         forAll(info, i)
01833         {
01834             if (info[i].hit())
01835             {
01836                 region[i] = triSurface::operator[](info[i].index()).region();
01837             }
01838             else
01839             {
01840                 region[i] = -1;
01841             }
01842         }
01843         return;
01844     }
01845 
01846 
01847     // Get query data (= local index of triangle)
01848     // ~~~~~~~~~~~~~~
01849 
01850     labelList triangleIndex(info.size());
01851     autoPtr<mapDistribute> mapPtr
01852     (
01853         calcLocalQueries
01854         (
01855             info,
01856             triangleIndex
01857         )
01858     );
01859     const mapDistribute& map = mapPtr();
01860 
01861 
01862     // Do my tests
01863     // ~~~~~~~~~~~
01864 
01865     const triSurface& s = static_cast<const triSurface&>(*this);
01866 
01867     region.setSize(triangleIndex.size());
01868 
01869     forAll(triangleIndex, i)
01870     {
01871         label triI = triangleIndex[i];
01872         region[i] = s[triI].region();
01873     }
01874 
01875 
01876     // Send back results
01877     // ~~~~~~~~~~~~~~~~~
01878 
01879     map.distribute
01880     (
01881         //Pstream::scheduled,
01882         //map.schedule            // note reverse schedule
01883         //(
01884         //    map.constructMap(),
01885         //    map.subMap()
01886         //),
01887         Pstream::nonBlocking,
01888         List<labelPair>(0),
01889         info.size(),
01890         map.constructMap(),     // what to send
01891         map.subMap(),           // what to receive
01892         region
01893     );
01894 }
01895 
01896 
01897 void Foam::distributedTriSurfaceMesh::getNormal
01898 (
01899     const List<pointIndexHit>& info,
01900     vectorField& normal
01901 ) const
01902 {
01903     if (!Pstream::parRun())
01904     {
01905         triSurfaceMesh::getNormal(info, normal);
01906         return;
01907     }
01908 
01909 
01910     // Get query data (= local index of triangle)
01911     // ~~~~~~~~~~~~~~
01912 
01913     labelList triangleIndex(info.size());
01914     autoPtr<mapDistribute> mapPtr
01915     (
01916         calcLocalQueries
01917         (
01918             info,
01919             triangleIndex
01920         )
01921     );
01922     const mapDistribute& map = mapPtr();
01923 
01924 
01925     // Do my tests
01926     // ~~~~~~~~~~~
01927 
01928     const triSurface& s = static_cast<const triSurface&>(*this);
01929 
01930     normal.setSize(triangleIndex.size());
01931 
01932     forAll(triangleIndex, i)
01933     {
01934         label triI = triangleIndex[i];
01935         normal[i] = s[triI].normal(s.points());
01936         normal[i] /= mag(normal[i]) + VSMALL;
01937     }
01938 
01939 
01940     // Send back results
01941     // ~~~~~~~~~~~~~~~~~
01942 
01943     map.distribute
01944     (
01945         //Pstream::scheduled,
01946         //map.schedule            // note reverse schedule
01947         //(
01948         //    map.constructMap(),
01949         //    map.subMap()
01950         //),
01951         Pstream::nonBlocking,
01952         List<labelPair>(0),
01953         info.size(),
01954         map.constructMap(),     // what to send
01955         map.subMap(),           // what to receive
01956         normal
01957     );
01958 }
01959 
01960 
01961 void Foam::distributedTriSurfaceMesh::getField
01962 (
01963     const List<pointIndexHit>& info,
01964     labelList& values
01965 ) const
01966 {
01967     if (!Pstream::parRun())
01968     {
01969         triSurfaceMesh::getField(info, values);
01970         return;
01971     }
01972 
01973     if (foundObject<triSurfaceLabelField>("values"))
01974     {
01975         const triSurfaceLabelField& fld = lookupObject<triSurfaceLabelField>
01976         (
01977             "values"
01978         );
01979 
01980 
01981         // Get query data (= local index of triangle)
01982         // ~~~~~~~~~~~~~~
01983 
01984         labelList triangleIndex(info.size());
01985         autoPtr<mapDistribute> mapPtr
01986         (
01987             calcLocalQueries
01988             (
01989                 info,
01990                 triangleIndex
01991             )
01992         );
01993         const mapDistribute& map = mapPtr();
01994 
01995 
01996         // Do my tests
01997         // ~~~~~~~~~~~
01998 
01999         values.setSize(triangleIndex.size());
02000 
02001         forAll(triangleIndex, i)
02002         {
02003             label triI = triangleIndex[i];
02004             values[i] = fld[triI];
02005         }
02006 
02007 
02008         // Send back results
02009         // ~~~~~~~~~~~~~~~~~
02010 
02011         map.distribute
02012         (
02013             Pstream::nonBlocking,
02014             List<labelPair>(0),
02015             info.size(),
02016             map.constructMap(),     // what to send
02017             map.subMap(),           // what to receive
02018             values
02019         );
02020     }
02021 }
02022 
02023 
02024 void Foam::distributedTriSurfaceMesh::getVolumeType
02025 (
02026     const pointField& points,
02027     List<volumeType>& volType
02028 ) const
02029 {
02030     FatalErrorIn
02031     (
02032         "distributedTriSurfaceMesh::getVolumeType"
02033         "(const pointField&, List<volumeType>&) const"
02034     )   << "Volume type not supported for distributed surfaces."
02035         << exit(FatalError);
02036 }
02037 
02038 
02039 // Subset the part of surface that is overlapping bb.
02040 Foam::triSurface Foam::distributedTriSurfaceMesh::overlappingSurface
02041 (
02042     const triSurface& s,
02043     const List<treeBoundBox>& bbs,
02044 
02045     labelList& subPointMap,
02046     labelList& subFaceMap
02047 )
02048 {
02049     // Determine what triangles to keep.
02050     boolList includedFace(s.size(), false);
02051 
02052     // Create a slightly larger bounding box.
02053     List<treeBoundBox> bbsX(bbs.size());
02054     const scalar eps = 1.0e-4;
02055     forAll(bbs, i)
02056     {
02057         const point mid = 0.5*(bbs[i].min() + bbs[i].max());
02058         const vector halfSpan = (1.0+eps)*(bbs[i].max() - mid);
02059 
02060         bbsX[i].min() = mid - halfSpan;
02061         bbsX[i].max() = mid + halfSpan;
02062     }
02063 
02064     forAll(s, triI)
02065     {
02066         const labelledTri& f = s[triI];
02067         const point& p0 = s.points()[f[0]];
02068         const point& p1 = s.points()[f[1]];
02069         const point& p2 = s.points()[f[2]];
02070 
02071         if (overlaps(bbsX, p0, p1, p2))
02072         {
02073             includedFace[triI] = true;
02074         }
02075     }
02076 
02077     return subsetMesh(s, includedFace, subPointMap, subFaceMap);
02078 }
02079 
02080 
02081 void Foam::distributedTriSurfaceMesh::distribute
02082 (
02083     const List<treeBoundBox>& bbs,
02084     const bool keepNonLocal,
02085     autoPtr<mapDistribute>& faceMap,
02086     autoPtr<mapDistribute>& pointMap
02087 )
02088 {
02089     // Get bbs of all domains
02090     // ~~~~~~~~~~~~~~~~~~~~~~
02091 
02092     {
02093         List<List<treeBoundBox> > newProcBb(Pstream::nProcs());
02094 
02095         switch(distType_)
02096         {
02097             case FOLLOW:
02098                 newProcBb[Pstream::myProcNo()].setSize(bbs.size());
02099                 forAll(bbs, i)
02100                 {
02101                     newProcBb[Pstream::myProcNo()][i] = bbs[i];
02102                 }
02103                 Pstream::gatherList(newProcBb);
02104                 Pstream::scatterList(newProcBb);
02105             break;
02106 
02107             case INDEPENDENT:
02108                 newProcBb = independentlyDistributedBbs(*this);
02109             break;
02110 
02111             case FROZEN:
02112                 return;
02113             break;
02114 
02115             default:
02116                 FatalErrorIn("distributedTriSurfaceMesh::distribute(..)")
02117                     << "Unsupported distribution type." << exit(FatalError);
02118             break;
02119         }
02120 
02121         //if (debug)
02122         //{
02123         //    Info<< "old bb:" << procBb_ << endl << endl;
02124         //    Info<< "new bb:" << newProcBb << endl << endl;
02125         //    Info<< "Same:" << (newProcBb == procBb_) << endl;
02126         //}
02127 
02128         if (newProcBb == procBb_)
02129         {
02130             return;
02131         }
02132         else
02133         {
02134             procBb_.transfer(newProcBb);
02135             dict_.set("bounds", procBb_[Pstream::myProcNo()]);
02136         }
02137     }
02138 
02139 
02140     // Debug information
02141     if (debug)
02142     {
02143         labelList nTris(Pstream::nProcs());
02144         nTris[Pstream::myProcNo()] = triSurface::size();
02145         Pstream::gatherList(nTris);
02146         Pstream::scatterList(nTris);
02147 
02148         Info<< "distributedTriSurfaceMesh::distribute : before distribution:"
02149             << endl
02150             << "\tproc\ttris" << endl;
02151 
02152         forAll(nTris, procI)
02153         {
02154             Info<< '\t' << procI << '\t' << nTris[procI] << endl;
02155         }
02156         Info<< endl;
02157     }
02158 
02159 
02160     // Use procBbs to determine which faces go where
02161     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
02162 
02163     labelListList faceSendMap(Pstream::nProcs());
02164     labelListList pointSendMap(Pstream::nProcs());
02165 
02166     forAll(procBb_, procI)
02167     {
02168         overlappingSurface
02169         (
02170             *this,
02171             procBb_[procI],
02172             pointSendMap[procI],
02173             faceSendMap[procI]
02174         );
02175 
02176         if (debug)
02177         {
02178             //Pout<< "Overlapping with proc " << procI
02179             //    << " faces:" << faceSendMap[procI].size()
02180             //    << " points:" << pointSendMap[procI].size() << endl << endl;
02181         }
02182     }
02183 
02184     if (keepNonLocal)
02185     {
02186         // Include in faceSendMap/pointSendMap the triangles that are
02187         // not mapped to any processor so they stay local.
02188 
02189         const triSurface& s = static_cast<const triSurface&>(*this);
02190 
02191         boolList includedFace(s.size(), true);
02192 
02193         forAll(faceSendMap, procI)
02194         {
02195             if (procI != Pstream::myProcNo())
02196             {
02197                 forAll(faceSendMap[procI], i)
02198                 {
02199                     includedFace[faceSendMap[procI][i]] = false;
02200                 }
02201             }
02202         }
02203 
02204         // Combine includedFace (all triangles that are not on any neighbour)
02205         // with overlap.
02206 
02207         forAll(faceSendMap[Pstream::myProcNo()], i)
02208         {
02209             includedFace[faceSendMap[Pstream::myProcNo()][i]] = true;
02210         }
02211 
02212         subsetMesh
02213         (
02214             s,
02215             includedFace,
02216             pointSendMap[Pstream::myProcNo()],
02217             faceSendMap[Pstream::myProcNo()]
02218         );
02219     }
02220 
02221 
02222     // Send over how many faces/points I need to receive
02223     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
02224 
02225     labelListList faceSendSizes(Pstream::nProcs());
02226     faceSendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
02227     forAll(faceSendMap, procI)
02228     {
02229         faceSendSizes[Pstream::myProcNo()][procI] = faceSendMap[procI].size();
02230     }
02231     Pstream::gatherList(faceSendSizes);
02232     Pstream::scatterList(faceSendSizes);
02233 
02234 
02235 
02236     // Exchange surfaces
02237     // ~~~~~~~~~~~~~~~~~
02238 
02239     // Storage for resulting surface
02240     List<labelledTri> allTris;
02241     pointField allPoints;
02242 
02243     labelListList faceConstructMap(Pstream::nProcs());
02244     labelListList pointConstructMap(Pstream::nProcs());
02245 
02246 
02247     // My own surface first
02248     // ~~~~~~~~~~~~~~~~~~~~
02249 
02250     {
02251         labelList pointMap;
02252         triSurface subSurface
02253         (
02254             subsetMesh
02255             (
02256                 *this,
02257                 faceSendMap[Pstream::myProcNo()],
02258                 pointMap
02259             )
02260         );
02261 
02262         allTris = subSurface;
02263         allPoints = subSurface.points();
02264 
02265         faceConstructMap[Pstream::myProcNo()] = identity
02266         (
02267             faceSendMap[Pstream::myProcNo()].size()
02268         );
02269         pointConstructMap[Pstream::myProcNo()] = identity
02270         (
02271             pointSendMap[Pstream::myProcNo()].size()
02272         );
02273     }
02274 
02275 
02276 
02277     // Send all
02278     // ~~~~~~~~
02279 
02280     forAll(faceSendSizes, procI)
02281     {
02282         if (procI != Pstream::myProcNo())
02283         {
02284             if (faceSendSizes[Pstream::myProcNo()][procI] > 0)
02285             {
02286                 OPstream str(Pstream::blocking, procI);
02287 
02288                 labelList pointMap;
02289                 triSurface subSurface
02290                 (
02291                     subsetMesh
02292                     (
02293                         *this,
02294                         faceSendMap[procI],
02295                         pointMap
02296                     )
02297                 );
02298 
02299                 //if (debug)
02300                 //{
02301                 //    Pout<< "Sending to " << procI
02302                 //        << " faces:" << faceSendMap[procI].size()
02303                 //        << " points:" << subSurface.points().size() << endl
02304                 //        << endl;
02305                 //}
02306 
02307                 str << subSurface;
02308            }
02309         }
02310     }
02311 
02312 
02313     // Receive and merge all
02314     // ~~~~~~~~~~~~~~~~~~~~~
02315 
02316     forAll(faceSendSizes, procI)
02317     {
02318         if (procI != Pstream::myProcNo())
02319         {
02320             if (faceSendSizes[procI][Pstream::myProcNo()] > 0)
02321             {
02322                 IPstream str(Pstream::blocking, procI);
02323 
02324                 // Receive
02325                 triSurface subSurface(str);
02326 
02327                 //if (debug)
02328                 //{
02329                 //    Pout<< "Received from " << procI
02330                 //        << " faces:" << subSurface.size()
02331                 //        << " points:" << subSurface.points().size() << endl
02332                 //        << endl;
02333                 //}
02334 
02335                 // Merge into allSurf
02336                 merge
02337                 (
02338                     mergeDist_,
02339                     subSurface,
02340                     subSurface.points(),
02341 
02342                     allTris,
02343                     allPoints,
02344                     faceConstructMap[procI],
02345                     pointConstructMap[procI]
02346                 );
02347 
02348                 //if (debug)
02349                 //{
02350                 //    Pout<< "Current merged surface : faces:" << allTris.size()
02351                 //        << " points:" << allPoints.size() << endl << endl;
02352                 //}
02353            }
02354         }
02355     }
02356 
02357 
02358     faceMap.reset
02359     (
02360         new mapDistribute
02361         (
02362             allTris.size(),
02363             faceSendMap,
02364             faceConstructMap,
02365             true
02366         )
02367     );
02368     pointMap.reset
02369     (
02370         new mapDistribute
02371         (
02372             allPoints.size(),
02373             pointSendMap,
02374             pointConstructMap,
02375             true
02376         )
02377     );
02378 
02379     // Construct triSurface. Reuse storage.
02380     triSurface::operator=(triSurface(allTris, patches(), allPoints, true));
02381 
02382     clearOut();
02383 
02384     // Regions stays same
02385     // Volume type stays same.
02386 
02387     distributeFields<label>(faceMap());
02388     distributeFields<scalar>(faceMap());
02389     distributeFields<vector>(faceMap());
02390     distributeFields<sphericalTensor>(faceMap());
02391     distributeFields<symmTensor>(faceMap());
02392     distributeFields<tensor>(faceMap());
02393 
02394     if (debug)
02395     {
02396         labelList nTris(Pstream::nProcs());
02397         nTris[Pstream::myProcNo()] = triSurface::size();
02398         Pstream::gatherList(nTris);
02399         Pstream::scatterList(nTris);
02400 
02401         Info<< "distributedTriSurfaceMesh::distribute : after distribution:"
02402             << endl
02403             << "\tproc\ttris" << endl;
02404 
02405         forAll(nTris, procI)
02406         {
02407             Info<< '\t' << procI << '\t' << nTris[procI] << endl;
02408         }
02409         Info<< endl;
02410     }
02411 }
02412 
02413 
02414 //- Write using given format, version and compression
02415 bool Foam::distributedTriSurfaceMesh::writeObject
02416 (
02417     IOstream::streamFormat fmt,
02418     IOstream::versionNumber ver,
02419     IOstream::compressionType cmp
02420 ) const
02421 {
02422     // Make sure dictionary goes to same directory as surface
02423     const_cast<fileName&>(dict_.instance()) = searchableSurface::instance();
02424 
02425     // Dictionary needs to be written in ascii - binary output not supported.
02426     return
02427         triSurfaceMesh::writeObject(fmt, ver, cmp)
02428      && dict_.writeObject(IOstream::ASCII, ver, cmp);
02429 }
02430 
02431 
02432 void Foam::distributedTriSurfaceMesh::writeStats(Ostream& os) const
02433 {
02434     boundBox bb;
02435     label nPoints;
02436     calcBounds(bb, nPoints);
02437     reduce(bb.min(), minOp<point>());
02438     reduce(bb.max(), maxOp<point>());
02439 
02440     os  << "Triangles    : " << returnReduce(triSurface::size(), sumOp<label>())
02441         << endl
02442         << "Vertices     : " << returnReduce(nPoints, sumOp<label>()) << endl
02443         << "Bounding Box : " << bb << endl;
02444 }
02445 
02446 
02447 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines