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

globalPoints.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 "globalPoints.H"
00027 #include <OpenFOAM/processorPolyPatch.H>
00028 #include <OpenFOAM/cyclicPolyPatch.H>
00029 #include <OpenFOAM/polyMesh.H>
00030 
00031 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
00032 
00033 defineTypeNameAndDebug(Foam::globalPoints, 0);
00034 
00035 
00036 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
00037 
00038 // Total number of points on processor patches. Is upper limit for number
00039 // of shared points
00040 Foam::label Foam::globalPoints::countPatchPoints
00041 (
00042     const polyBoundaryMesh& patches
00043 )
00044 {
00045     label nTotPoints = 0;
00046 
00047     forAll(patches, patchI)
00048     {
00049         const polyPatch& pp = patches[patchI];
00050 
00051         if
00052         (
00053             (Pstream::parRun() && isA<processorPolyPatch>(pp))
00054          || isA<cyclicPolyPatch>(pp)
00055         )
00056         {
00057             nTotPoints += pp.nPoints();
00058         }
00059     }
00060 
00061     return nTotPoints;
00062 }
00063 
00064 
00065 // Collect all topological information about a point on a patch.
00066 // (this information is the patch faces using the point and the relative
00067 // position of the point in the face)
00068 void Foam::globalPoints::addToSend
00069 (
00070     const primitivePatch& pp,
00071     const label patchPointI,
00072     const procPointList& knownInfo,
00073 
00074     DynamicList<label>& patchFaces,
00075     DynamicList<label>& indexInFace,
00076     DynamicList<procPointList>& allInfo
00077 )
00078 {
00079     label meshPointI = pp.meshPoints()[patchPointI];
00080 
00081     // Add all faces using the point so we are sure we find it on the
00082     // other side.
00083     const labelList& pFaces = pp.pointFaces()[patchPointI];
00084 
00085     forAll(pFaces, i)
00086     {
00087         label patchFaceI = pFaces[i];
00088 
00089         const face& f = pp[patchFaceI];
00090 
00091         patchFaces.append(patchFaceI);
00092         indexInFace.append(findIndex(f, meshPointI));
00093         allInfo.append(knownInfo);
00094     }
00095 }
00096 
00097 
00098 // Add nbrInfo to myInfo. Return true if anything changed.
00099 // nbrInfo is for a point a list of all the processors using it (and the
00100 // meshPoint label on that processor)
00101 bool Foam::globalPoints::mergeInfo
00102 (
00103     const procPointList& nbrInfo,
00104     procPointList& myInfo
00105 )
00106 {
00107     // Indices of entries in nbrInfo not yet in myInfo.
00108     DynamicList<label> newInfo(nbrInfo.size());
00109 
00110     forAll(nbrInfo, i)
00111     {
00112         const procPoint& info = nbrInfo[i];
00113 
00114         // Check if info already in myInfo.
00115         label index = -1;
00116 
00117         forAll(myInfo, j)
00118         {
00119             if (myInfo[j] == info)
00120             {
00121                 // Already have information for processor/point combination
00122                 // in my list so skip.
00123                 index = j;
00124 
00125                 break;
00126             }
00127         }
00128 
00129         if (index == -1)
00130         {
00131             // Mark this information as being not yet in myInfo
00132             newInfo.append(i);
00133         }
00134     }
00135 
00136     newInfo.shrink();
00137 
00138     // Append all nbrInfos referenced in newInfo to myInfo.
00139 
00140     label index = myInfo.size();
00141 
00142     myInfo.setSize(index + newInfo.size());
00143 
00144     forAll(newInfo, i)
00145     {
00146         myInfo[index++] = nbrInfo[newInfo[i]];
00147     }
00148 
00149     // Did anything change?
00150     return newInfo.size() > 0;
00151 }
00152 
00153 
00154 // Updates database of current information on meshpoints with nbrInfo.
00155 // Uses mergeInfo above. Returns true if data kept for meshPointI changed.
00156 bool Foam::globalPoints::storeInfo
00157 (
00158     const procPointList& nbrInfo,
00159     const label meshPointI
00160 )
00161 {
00162     label infoChanged = false;
00163 
00164     // Get the index into the procPoints list.
00165     Map<label>::iterator iter = meshToProcPoint_.find(meshPointI);
00166 
00167     if (iter != meshToProcPoint_.end())
00168     {
00169         procPointList& knownInfo = procPoints_[iter()];
00170 
00171         if (mergeInfo(nbrInfo, knownInfo))
00172         {
00173             infoChanged = true;
00174         }
00175     }
00176     else
00177     {
00178         procPointList knownInfo(1);
00179         knownInfo[0][0] = Pstream::myProcNo();
00180         knownInfo[0][1] = meshPointI;
00181 
00182         if (mergeInfo(nbrInfo, knownInfo))
00183         {
00184             // Update addressing from into procPoints
00185             meshToProcPoint_.insert(meshPointI, procPoints_.size());
00186             // Insert into list of equivalences.
00187             procPoints_.append(knownInfo);
00188 
00189             infoChanged = true;
00190         }
00191     }
00192     return infoChanged;
00193 }
00194 
00195 
00196 // Insert my own points into structure and mark as changed.
00197 void Foam::globalPoints::initOwnPoints
00198 (
00199     const bool allPoints,
00200     labelHashSet& changedPoints
00201 )
00202 {
00203     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
00204 
00205     forAll(patches, patchI)
00206     {
00207         const polyPatch& pp = patches[patchI];
00208 
00209         if
00210         (
00211             (Pstream::parRun() && isA<processorPolyPatch>(pp))
00212          || isA<cyclicPolyPatch>(pp)
00213         )
00214         {
00215             const labelList& meshPoints = pp.meshPoints();
00216 
00217             if (allPoints)
00218             {
00219                 // All points on patch
00220                 forAll(meshPoints, i)
00221                 {
00222                     label meshPointI = meshPoints[i];
00223 
00224                     procPointList knownInfo(1);
00225                     knownInfo[0][0] = Pstream::myProcNo();
00226                     knownInfo[0][1] = meshPointI;
00227 
00228                     // Update addressing from meshpoint to index in procPoints
00229                     meshToProcPoint_.insert(meshPointI, procPoints_.size());
00230                     // Insert into list of equivalences.
00231                     procPoints_.append(knownInfo);
00232 
00233                     // Update changedpoints info.
00234                     changedPoints.insert(meshPointI);
00235                 }
00236             }
00237             else
00238             {
00239                 // Boundary points only
00240                 const labelList& boundaryPoints = pp.boundaryPoints();
00241 
00242                 forAll(boundaryPoints, i)
00243                 {
00244                     label meshPointI = meshPoints[boundaryPoints[i]];
00245 
00246                     procPointList knownInfo(1);
00247                     knownInfo[0][0] = Pstream::myProcNo();
00248                     knownInfo[0][1] = meshPointI;
00249 
00250                     // Update addressing from meshpoint to index in procPoints
00251                     meshToProcPoint_.insert(meshPointI, procPoints_.size());
00252                     // Insert into list of equivalences.
00253                     procPoints_.append(knownInfo);
00254 
00255                     // Update changedpoints info.
00256                     changedPoints.insert(meshPointI);
00257                 }
00258             }
00259         }
00260     }
00261 }
00262 
00263 
00264 // Send all my info on changedPoints_ to my neighbours.
00265 void Foam::globalPoints::sendPatchPoints(const labelHashSet& changedPoints)
00266  const
00267 {
00268     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
00269 
00270     forAll(patches, patchI)
00271     {
00272         const polyPatch& pp = patches[patchI];
00273 
00274         if (Pstream::parRun() && isA<processorPolyPatch>(pp))
00275         {
00276             // Information to send:
00277             // patch face
00278             DynamicList<label> patchFaces(pp.nPoints());
00279             // index in patch face
00280             DynamicList<label> indexInFace(pp.nPoints());
00281             // all information I currently hold about this patchPoint
00282             DynamicList<procPointList> allInfo(pp.nPoints());
00283 
00284 
00285             // Now collect information on all mesh points mentioned in
00286             // changedPoints. Note that these points only should occur on
00287             // processorPatches (or rather this is a limitation!).
00288 
00289             const labelList& meshPoints = pp.meshPoints();
00290 
00291             forAll(meshPoints, patchPointI)
00292             {
00293                 label meshPointI = meshPoints[patchPointI];
00294 
00295                 if (changedPoints.found(meshPointI))
00296                 {
00297                     label index = meshToProcPoint_[meshPointI];
00298 
00299                     const procPointList& knownInfo = procPoints_[index];
00300 
00301                     // Add my information about meshPointI to the send buffers
00302                     addToSend
00303                     (
00304                         pp,
00305                         patchPointI,
00306                         knownInfo,
00307 
00308                         patchFaces,
00309                         indexInFace,
00310                         allInfo
00311                     );
00312                 }
00313             }
00314             patchFaces.shrink();
00315             indexInFace.shrink();
00316             allInfo.shrink();
00317 
00318             // Send to neighbour
00319             {
00320                 const processorPolyPatch& procPatch =
00321                     refCast<const processorPolyPatch>(pp);
00322 
00323                 if (debug)
00324                 {
00325                     Pout<< " Sending to "
00326                         << procPatch.neighbProcNo() << "   point information:"
00327                         << patchFaces.size() << endl;
00328                 }
00329 
00330                 OPstream toNeighbour
00331                 (
00332                     Pstream::blocking,
00333                     procPatch.neighbProcNo()
00334                 );
00335 
00336                 toNeighbour << patchFaces << indexInFace << allInfo;
00337             }
00338         }
00339     }
00340 }
00341 
00342 
00343 // Receive all my neighbours' information and merge with mine.
00344 // After finishing will have updated
00345 // - procPoints_ : all neighbour information merged in.
00346 // - meshToProcPoint_
00347 // - changedPoints: all meshPoints for which something changed.
00348 void Foam::globalPoints::receivePatchPoints(labelHashSet& changedPoints)
00349 {
00350     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
00351 
00352     // Reset changed points
00353     changedPoints.clear();
00354 
00355     forAll(patches, patchI)
00356     {
00357         const polyPatch& pp = patches[patchI];
00358 
00359         if (Pstream::parRun() && isA<processorPolyPatch>(pp))
00360         {
00361             const processorPolyPatch& procPatch =
00362                 refCast<const processorPolyPatch>(pp);
00363 
00364             labelList patchFaces;
00365             labelList indexInFace;
00366             List<procPointList> nbrInfo;
00367 
00368             {
00369                 IPstream fromNeighbour
00370                 (
00371                     Pstream::blocking,
00372                     procPatch.neighbProcNo()
00373                 );
00374                 fromNeighbour >> patchFaces >> indexInFace >> nbrInfo;
00375             }
00376 
00377             if (debug)
00378             {
00379                 Pout<< " Received from "
00380                     << procPatch.neighbProcNo() << "   point information:"
00381                     << patchFaces.size() << endl;
00382             }
00383 
00384             forAll(patchFaces, i)
00385             {
00386                 const face& f = pp[patchFaces[i]];
00387 
00388                 // Get index in this face from index on face on other side.
00389                 label index = (f.size() - indexInFace[i]) % f.size();
00390 
00391                 // Get the meshpoint on my side
00392                 label meshPointI = f[index];
00393 
00394                 //const procPointList& info = nbrInfo[i];
00395                 //Pout << "Received for my coord "
00396                 //    << mesh_.points()[meshPointI];
00397                 //
00398                 //forAll(info, j)
00399                 //{
00400                 //    Pout<< ' ' <<info[j];
00401                 //}
00402                 //Pout<< endl;
00403 
00404                 if (storeInfo(nbrInfo[i], meshPointI))
00405                 {
00406                     changedPoints.insert(meshPointI);
00407                 }
00408             }
00409         }
00410         else if (isA<cyclicPolyPatch>(pp))
00411         {
00412             // Handle cyclics: send lower half to upper half and vice versa.
00413             // Or since they both are in memory just do it point by point.
00414 
00415             const cyclicPolyPatch& cycPatch =
00416                 refCast<const cyclicPolyPatch>(pp);
00417 
00418             const labelList& meshPoints = pp.meshPoints();
00419 
00420             //const edgeList& connections = cycPatch.coupledPoints();
00421             const edgeList connections(coupledPoints(cycPatch));
00422 
00423             forAll(connections, i)
00424             {
00425                 const edge& e = connections[i];
00426 
00427                 label meshPointA = meshPoints[e[0]];
00428                 label meshPointB = meshPoints[e[1]];
00429 
00430                 // Do we have information on meshPointA?
00431                 Map<label>::iterator procPointA =
00432                     meshToProcPoint_.find(meshPointA);
00433 
00434                 if (procPointA != meshToProcPoint_.end())
00435                 {
00436                     // Store A info onto pointB
00437                     if (storeInfo(procPoints_[procPointA()], meshPointB))
00438                     {
00439                         changedPoints.insert(meshPointB);
00440                     }
00441                 }
00442 
00443                 // Same for info on pointB
00444                 Map<label>::iterator procPointB =
00445                     meshToProcPoint_.find(meshPointB);
00446 
00447                 if (procPointB != meshToProcPoint_.end())
00448                 {
00449                     // Store B info onto pointA
00450                     if (storeInfo(procPoints_[procPointB()], meshPointA))
00451                     {
00452                         changedPoints.insert(meshPointA);
00453                     }
00454                 }
00455             }
00456         }
00457     }
00458 }
00459 
00460 
00461 // Remove entries which are handled by normal face-face communication. I.e.
00462 // those points where the equivalence list is only me and my (face)neighbour
00463 void Foam::globalPoints::remove(const Map<label>& directNeighbours)
00464 {
00465     // Save old ones.
00466     Map<label> oldMeshToProcPoint(meshToProcPoint_);
00467     meshToProcPoint_.clear();
00468 
00469     List<procPointList> oldProcPoints;
00470     oldProcPoints.transfer(procPoints_);
00471 
00472     // Go through all equivalences
00473     for
00474     (
00475         Map<label>::const_iterator iter = oldMeshToProcPoint.begin();
00476         iter != oldMeshToProcPoint.end();
00477         ++iter
00478     )
00479     {
00480         label meshPointI = iter.key();
00481         const procPointList& pointInfo = oldProcPoints[iter()];
00482 
00483         if (pointInfo.size() == 2)
00484         {
00485             // I will be in this equivalence list.
00486             // Check whether my direct (=face) neighbour
00487             // is in it. This would be an ordinary connection and can be
00488             // handled by normal face-face connectivity.
00489 
00490             const procPoint& a = pointInfo[0];
00491             const procPoint& b = pointInfo[1];
00492 
00493             if
00494             (
00495                 (a[0] == Pstream::myProcNo() && directNeighbours.found(a[1]))
00496              || (b[0] == Pstream::myProcNo() && directNeighbours.found(b[1]))
00497             )
00498             {
00499                 // Normal faceNeighbours
00500                 if (a[0] == Pstream::myProcNo())
00501                 {
00502                     //Pout<< "Removing direct neighbour:"
00503                     //    << mesh_.points()[a[1]]
00504                     //    << endl;
00505                 }
00506                 else if (b[0] == Pstream::myProcNo())
00507                 {
00508                     //Pout<< "Removing direct neighbour:"
00509                     //    << mesh_.points()[b[1]]
00510                     //    << endl;
00511                 }
00512             }
00513             else
00514             {
00515                 // This condition will be very rare: points are used by
00516                 // two processors which are not face-face connected.
00517                 // e.g.
00518                 // +------+------+
00519                 // | wall |  B   |
00520                 // +------+------+
00521                 // |   A  | wall |
00522                 // +------+------+
00523                 // Processor A and B share a point. Note that this only will
00524                 // be found if the two domains are face connected at all
00525                 // (not shown in the picture)
00526 
00527                 meshToProcPoint_.insert(meshPointI, procPoints_.size());
00528                 procPoints_.append(pointInfo);
00529             }
00530         }
00531         else if (pointInfo.size() == 1)
00532         {
00533             // This happens for 'wedge' like cyclics where the two halves
00534             // come together in the same point so share the same meshPoint.
00535             // So this meshPoint will have info of size one only.
00536             if
00537             (
00538                 pointInfo[0][0] != Pstream::myProcNo()
00539              || !directNeighbours.found(pointInfo[0][1])
00540             )
00541             {
00542                 meshToProcPoint_.insert(meshPointI, procPoints_.size());
00543                 procPoints_.append(pointInfo);
00544             }
00545         }
00546         else
00547         {
00548             meshToProcPoint_.insert(meshPointI, procPoints_.size());
00549             procPoints_.append(pointInfo);
00550         }
00551     }
00552 
00553     procPoints_.shrink();
00554 }
00555 
00556 
00557 // Get (indices of) points for which I am master (= lowest numbered point on
00558 // lowest numbered processor).
00559 // (equivalence lists should be complete by now)
00560 Foam::labelList Foam::globalPoints::getMasterPoints() const
00561 {
00562     labelList masterPoints(nPatchPoints_);
00563     label nMaster = 0;
00564 
00565     // Go through all equivalences and determine meshPoints where I am master.
00566     for
00567     (
00568         Map<label>::const_iterator iter = meshToProcPoint_.begin();
00569         iter != meshToProcPoint_.end();
00570         ++iter
00571     )
00572     {
00573         label meshPointI = iter.key();
00574         const procPointList& pointInfo = procPoints_[iter()];
00575 
00576         if (pointInfo.size() < 2)
00577         {
00578             // Points should have an equivalence list >= 2 since otherwise
00579             // they would be direct neighbours and have been removed in the
00580             // call to 'remove'.
00581             Pout<< "MeshPoint:" << meshPointI
00582                 << " coord:" << mesh_.points()[meshPointI]
00583                 << " has no corresponding point on a neighbouring processor"
00584                 << endl;
00585             FatalErrorIn("globalPoints::getMasterPoints()")
00586                 << '[' << Pstream::myProcNo() << ']'
00587                 << " MeshPoint:" << meshPointI
00588                 << " coord:" << mesh_.points()[meshPointI]
00589                 << " has no corresponding point on a neighbouring processor"
00590                 << abort(FatalError);
00591         }
00592         else
00593         {
00594             // Find lowest processor and lowest mesh point on this processor.
00595             label lowestProcI = labelMax;
00596             label lowestPointI = labelMax;
00597 
00598             forAll(pointInfo, i)
00599             {
00600                 label proc = pointInfo[i][0];
00601 
00602                 if (proc < lowestProcI)
00603                 {
00604                     lowestProcI = proc;
00605                     lowestPointI = pointInfo[i][1];
00606                 }
00607                 else if (proc == lowestProcI)
00608                 {
00609                     lowestPointI = min(lowestPointI, pointInfo[i][1]);
00610                 }
00611             }
00612 
00613             if
00614             (
00615                 lowestProcI == Pstream::myProcNo()
00616              && lowestPointI == meshPointI
00617             )
00618             {
00619                 // I am lowest numbered processor and point. Add to my list.
00620                 masterPoints[nMaster++] = meshPointI;
00621             }
00622         }
00623     }
00624 
00625     masterPoints.setSize(nMaster);
00626 
00627     return masterPoints;
00628 }
00629 
00630 
00631 // Send subset of lists
00632 void Foam::globalPoints::sendSharedPoints(const labelList& changedIndices) const
00633 {
00634     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
00635 
00636     forAll(patches, patchI)
00637     {
00638         const polyPatch& pp = patches[patchI];
00639 
00640         if (Pstream::parRun() && isA<processorPolyPatch>(pp))
00641         {
00642             const processorPolyPatch& procPatch =
00643                 refCast<const processorPolyPatch>(pp);
00644 
00645             OPstream toNeighbour(Pstream::blocking, procPatch.neighbProcNo());
00646 
00647             if (debug)
00648             {
00649                 Pout<< "Sending to " << procPatch.neighbProcNo()
00650                     << "  changed sharedPoints info:"
00651                     << changedIndices.size() << endl;
00652             }
00653 
00654             toNeighbour
00655                 << UIndirectList<label>(sharedPointAddr_, changedIndices)()
00656                 << UIndirectList<label>(sharedPointLabels_, changedIndices)();
00657         }
00658     }
00659 }
00660 
00661 
00662 // Receive shared point indices for all my shared points. Note that since
00663 // there are only a few here we can build a reverse map using the meshpoint
00664 // instead of doing all this relative point indexing (patch face + index in
00665 // face) as in send/receivePatchPoints
00666 void Foam::globalPoints::receiveSharedPoints(labelList& changedIndices)
00667 {
00668     changedIndices.setSize(sharedPointAddr_.size());
00669     label nChanged = 0;
00670 
00671     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
00672 
00673     // Receive and set shared points
00674     forAll(patches, patchI)
00675     {
00676         const polyPatch& pp = patches[patchI];
00677 
00678         if (Pstream::parRun() && isA<processorPolyPatch>(pp))
00679         {
00680             const processorPolyPatch& procPatch =
00681                 refCast<const processorPolyPatch>(pp);
00682 
00683             // Map from neighbouring meshPoint to sharedPoint)
00684             Map<label> nbrSharedPoints(sharedPointAddr_.size());
00685 
00686             {
00687                 // Receive meshPoints on neighbour and sharedPoints and build
00688                 // map from it. Note that we could have built the map on the
00689                 // neighbour and sent it over.
00690                 labelList nbrSharedPointAddr;
00691                 labelList nbrSharedPointLabels;
00692 
00693                 {
00694                     IPstream fromNeighbour
00695                     (
00696                         Pstream::blocking,
00697                         procPatch.neighbProcNo()
00698                     );
00699                     fromNeighbour >> nbrSharedPointAddr >> nbrSharedPointLabels;
00700                 }
00701 
00702                 // Insert into to map
00703                 forAll(nbrSharedPointLabels, i)
00704                 {
00705                     nbrSharedPoints.insert
00706                     (
00707                         nbrSharedPointLabels[i], // meshpoint on neighbour
00708                         nbrSharedPointAddr[i]    // sharedPoint label
00709                     );
00710                 }
00711             }
00712 
00713 
00714             // Merge into whatever information I hold.
00715             for
00716             (
00717                 Map<label>::const_iterator iter = meshToProcPoint_.begin();
00718                 iter != meshToProcPoint_.end();
00719                 ++iter
00720             )
00721             {
00722                 label meshPointI = iter.key();
00723                 label index = iter();
00724 
00725                 if (sharedPointAddr_[index] == -1)
00726                 {
00727                     // No shared point known yet for this meshPoint.
00728                     // See if was received from neighbour.
00729                     const procPointList& knownInfo = procPoints_[index];
00730 
00731                     // Check through the whole equivalence list for any
00732                     // point from the neighbour.
00733                     forAll(knownInfo, j)
00734                     {
00735                         const procPoint& info = knownInfo[j];
00736 
00737                         if
00738                         (
00739                             (info[0] == procPatch.neighbProcNo())
00740                          && nbrSharedPoints.found(info[1])
00741                         )
00742                         {
00743                             // So this knownInfo contains the neighbour point
00744                             label sharedPointI = nbrSharedPoints[info[1]];
00745 
00746                             sharedPointAddr_[index] = sharedPointI;
00747                             sharedPointLabels_[index] = meshPointI;
00748                             changedIndices[nChanged++] = index;
00749 
00750                             break;
00751                         }
00752                     }
00753                 }
00754             }
00755         }
00756         else if (isA<cyclicPolyPatch>(pp))
00757         {
00758             const cyclicPolyPatch& cycPatch =
00759                 refCast<const cyclicPolyPatch>(pp);
00760 
00761             // Build map from meshPoint to sharedPoint
00762             Map<label> meshToSharedPoint(sharedPointAddr_.size());
00763             forAll(sharedPointLabels_, i)
00764             {
00765                 label meshPointI = sharedPointLabels_[i];
00766 
00767                 meshToSharedPoint.insert(meshPointI, sharedPointAddr_[i]);
00768             }
00769 
00770             // Sync all info.
00771             //const edgeList& connections = cycPatch.coupledPoints();
00772             const edgeList connections(coupledPoints(cycPatch));
00773 
00774             forAll(connections, i)
00775             {
00776                 const edge& e = connections[i];
00777 
00778                 label meshPointA = pp.meshPoints()[e[0]];
00779                 label meshPointB = pp.meshPoints()[e[1]];
00780 
00781                 // Do we already have shared point for meshPointA?
00782                 Map<label>::iterator fndA = meshToSharedPoint.find(meshPointA);
00783                 Map<label>::iterator fndB = meshToSharedPoint.find(meshPointB);
00784 
00785                 if (fndA != meshToSharedPoint.end())
00786                 {
00787                     if (fndB != meshToSharedPoint.end())
00788                     {
00789                         if (fndA() != fndB())
00790                         {
00791                             FatalErrorIn
00792                             (
00793                                 "globalPoints::receiveSharedPoints"
00794                                 "(labelList&)"
00795                             )   << "On patch " << pp.name()
00796                                 << " connected points " << meshPointA
00797                                 << ' ' << mesh_.points()[meshPointA]
00798                                 << " and " << meshPointB
00799                                 << ' ' << mesh_.points()[meshPointB]
00800                                 << " are mapped to different shared points: "
00801                                 << fndA() << " and " << fndB()
00802                                 << abort(FatalError);
00803                         }
00804                     }
00805                     else
00806                     {
00807                         // No shared point yet for B.
00808                         label sharedPointI = fndA();
00809 
00810                         // Store shared point for meshPointB
00811                         label index = meshToProcPoint_[meshPointB];
00812 
00813                         sharedPointAddr_[index] = sharedPointI;
00814                         sharedPointLabels_[index] = meshPointB;
00815                         changedIndices[nChanged++] = index;
00816                     }
00817                 }
00818                 else
00819                 {
00820                     // No shared point yet for A.
00821                     if (fndB != meshToSharedPoint.end())
00822                     {
00823                         label sharedPointI = fndB();
00824 
00825                         // Store shared point for meshPointA
00826                         label index = meshToProcPoint_[meshPointA];
00827 
00828                         sharedPointAddr_[index] = sharedPointI;
00829                         sharedPointLabels_[index] = meshPointA;
00830                         changedIndices[nChanged++] = index;
00831                     }
00832                 }
00833             }
00834         }
00835     }
00836 
00837     changedIndices.setSize(nChanged);
00838 }
00839 
00840 
00841 Foam::edgeList Foam::globalPoints::coupledPoints(const cyclicPolyPatch& pp)
00842 {
00843     // Look at cyclic patch as two halves, A and B.
00844     // Now all we know is that relative face index in halfA is same
00845     // as coupled face in halfB and also that the 0th vertex
00846     // corresponds.
00847 
00848     // From halfA point to halfB or -1.
00849     labelList coupledPoint(pp.nPoints(), -1);
00850 
00851     for (label patchFaceA = 0; patchFaceA < pp.size()/2; patchFaceA++)
00852     {
00853         const face& fA = pp.localFaces()[patchFaceA];
00854 
00855         forAll(fA, indexA)
00856         {
00857             label patchPointA = fA[indexA];
00858 
00859             if (coupledPoint[patchPointA] == -1)
00860             {
00861                 const face& fB = pp.localFaces()[patchFaceA + pp.size()/2];
00862 
00863                 label indexB = (fB.size() - indexA) % fB.size();
00864 
00865                 coupledPoint[patchPointA] = fB[indexB];
00866             }
00867         }
00868     }
00869 
00870     edgeList connected(pp.nPoints());
00871 
00872     // Extract coupled points.
00873     label connectedI = 0;
00874 
00875     forAll(coupledPoint, i)
00876     {
00877         if (coupledPoint[i] != -1)
00878         {
00879             connected[connectedI++] = edge(i, coupledPoint[i]);
00880         }
00881     }
00882 
00883     connected.setSize(connectedI);
00884 
00885     return connected;
00886 }
00887 
00888 
00889 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
00890 
00891 // Construct from components
00892 Foam::globalPoints::globalPoints(const polyMesh& mesh)
00893 :
00894     mesh_(mesh),
00895     nPatchPoints_(countPatchPoints(mesh.boundaryMesh())),
00896     procPoints_(nPatchPoints_),
00897     meshToProcPoint_(nPatchPoints_),
00898     sharedPointAddr_(0),
00899     sharedPointLabels_(0),
00900     nGlobalPoints_(0)
00901 {
00902     if (debug)
00903     {
00904         Pout<< "globalPoints::globalPoints(const polyMesh&) : "
00905             << "doing processor to processor communication to get sharedPoints"
00906             << endl;
00907     }
00908 
00909     labelHashSet changedPoints(nPatchPoints_);
00910 
00911     // Initialize procPoints with my patch points. Keep track of points
00912     // inserted (in changedPoints)
00913     // There are two possible forms of this:
00914     // - initialize with all patch points (allPoints = true). This causes all
00915     //   patch points to be exchanged so a lot of information gets stored and
00916     //   transferred. This all gets filtered out later when removing the
00917     //   equivalence lists of size 2.
00918     // - initialize with boundary points of patches only (allPoints = false).
00919     //   This should work for all decompositions except extreme ones where a
00920     //   shared point is not on the boundary of any processor patches using it.
00921     //   This would happen if a domain was pinched such that two patches share
00922     //   a point or edge.
00923     initOwnPoints(true, changedPoints);
00924 
00925     // Do one exchange iteration to get neighbour points.
00926     sendPatchPoints(changedPoints);
00927     receivePatchPoints(changedPoints);
00928 
00929 
00930     // Save neighbours reachable through face-face communication.
00931     Map<label> neighbourList(meshToProcPoint_);
00932 
00933 
00934     // Exchange until nothing changes on all processors.
00935     bool changed = false;
00936 
00937     do
00938     {
00939         sendPatchPoints(changedPoints);
00940         receivePatchPoints(changedPoints);
00941 
00942         changed = changedPoints.size() > 0;
00943         reduce(changed, orOp<bool>());
00944 
00945     } while (changed);
00946 
00947 
00948     // Remove direct neighbours from point equivalences.
00949     remove(neighbourList);
00950 
00951 
00952     //Pout<< "After removing locally connected points:" << endl;
00953     //for
00954     //(
00955     //    Map<label>::const_iterator iter = meshToProcPoint_.begin();
00956     //    iter != meshToProcPoint_.end();
00957     //    ++iter
00958     //)
00959     //{
00960     //    label meshPointI = iter.key();
00961     //    const procPointList& pointInfo = procPoints_[iter()];
00962     //
00963     //    forAll(pointInfo, i)
00964     //    {
00965     //        Pout<< "    pointI:" << meshPointI << ' '
00966     //            << mesh.points()[meshPointI]
00967     //            << " connected to proc " << pointInfo[i][0]
00968     //            << " point:" << pointInfo[i][1]
00969     //        << endl;
00970     //    }
00971     //}
00972 
00973 
00974     // We now have - in procPoints_ - a list of points which are shared between
00975     // multiple processors. These are the ones for which are sharedPoint
00976     // needs to be determined. This is done by having the lowest numbered
00977     // processor in the equivalence list 'ask' for a sharedPoint number
00978     // and then distribute it across processor patches to the non-master
00979     // processors. Note: below piece of coding is not very efficient. Uses
00980     // a Map where possibly it shouldn't
00981 
00982     // Initialize sharedPoint addressing. Is for every entry in procPoints_
00983     // the sharedPoint.
00984     sharedPointAddr_.setSize(meshToProcPoint_.size());
00985     sharedPointAddr_ = -1;
00986     sharedPointLabels_.setSize(meshToProcPoint_.size());
00987     sharedPointLabels_ = -1;
00988 
00989 
00990     // Get points for which I am master (lowest numbered proc)
00991     labelList masterPoints(getMasterPoints());
00992 
00993 
00994     // Determine number of master points on all processors.
00995     labelList sharedPointSizes(Pstream::nProcs());
00996     sharedPointSizes[Pstream::myProcNo()] = masterPoints.size();
00997 
00998     Pstream::gatherList(sharedPointSizes);
00999     Pstream::scatterList(sharedPointSizes);
01000 
01001     if (debug)
01002     {
01003         Pout<< "sharedPointSizes:" << sharedPointSizes << endl;
01004     }
01005 
01006     // Get total number of shared points
01007     nGlobalPoints_ = 0;
01008     forAll(sharedPointSizes, procI)
01009     {
01010         nGlobalPoints_ += sharedPointSizes[procI];
01011     }
01012 
01013     // Assign sharedPoint labels. Every processor gets assigned consecutive
01014     // numbers for its master points.
01015     // These are assigned in processor order so processor0 gets
01016     // 0..masterPoints.size()-1 etc.
01017 
01018     // My point labels start after those of lower numbered processors
01019     label sharedPointI = 0;
01020     for (label procI = 0; procI < Pstream::myProcNo(); procI++)
01021     {
01022         sharedPointI += sharedPointSizes[procI];
01023     }
01024 
01025     forAll(masterPoints, i)
01026     {
01027         label meshPointI = masterPoints[i];
01028 
01029         label index = meshToProcPoint_[meshPointI];
01030 
01031         sharedPointLabels_[index] = meshPointI;
01032         sharedPointAddr_[index] = sharedPointI++;
01033     }
01034 
01035 
01036     // Now we have a sharedPointLabel for some of the entries in procPoints.
01037     // Send this information to neighbours. Receive their information.
01038     // Loop until nothing changes.
01039 
01040     // Initial subset to send is points for which I have sharedPoints
01041     labelList changedIndices(sharedPointAddr_.size());
01042     label nChanged = 0;
01043 
01044     forAll(sharedPointAddr_, i)
01045     {
01046         if (sharedPointAddr_[i] != -1)
01047         {
01048             changedIndices[nChanged++] = i;
01049         }
01050     }
01051     changedIndices.setSize(nChanged);
01052 
01053     changed = false;
01054 
01055     do
01056     {
01057         if (debug)
01058         {
01059             Pout<< "Determined " << changedIndices.size() << " shared points."
01060                 << " Exchanging them" << endl;
01061         }
01062         sendSharedPoints(changedIndices);
01063         receiveSharedPoints(changedIndices);
01064 
01065         changed = changedIndices.size() > 0;
01066         reduce(changed, orOp<bool>());
01067 
01068     } while (changed);
01069 
01070 
01071     forAll(sharedPointLabels_, i)
01072     {
01073         if (sharedPointLabels_[i] == -1)
01074         {
01075             FatalErrorIn("globalPoints::globalPoints(const polyMesh& mesh)")
01076                 << "Problem: shared point on processor " << Pstream::myProcNo()
01077                 << " not set at index " << sharedPointLabels_[i] << endl
01078                 << "This might mean the individual processor domains are not"
01079                 << " connected and the overall domain consists of multiple"
01080                 << " regions. You can check this with checkMesh"
01081                 << abort(FatalError);
01082         }
01083     }
01084 
01085     if (debug)
01086     {
01087         Pout<< "globalPoints::globalPoints(const polyMesh&) : "
01088             << "Finished global points" << endl;
01089     }
01090 }
01091 
01092 
01093 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines