src/mesh/cfMesh/meshLibrary/utilities/helperFunctions/helperFunctionsFrontalMarking.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | cfMesh: A library for mesh generation
   4    \\    /   O peration     |
   5     \\  /    A nd           | Author: Franjo Juretic (franjo.juretic@c-fields.com)
   6      \\/     M anipulation  | Copyright (C) Creative Fields, Ltd.
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of cfMesh.
  10
  11     cfMesh is free software; you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation; either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     cfMesh is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with cfMesh.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Description
  25
  26 \*---------------------------------------------------------------------------*/
  27
  28 #include "error.H"
  29 #include "helperFunctionsFrontalMarking.H"
  30 #include "DynList.H"
  31 #include "labelPair.H"
  32 #include "HashSet.H"
  33
  34 # ifdef USE_OMP
  35 #include <omp.h>
  36 # endif
  37
  38 namespace Foam
  39 {
  40
  41 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *//
  42
  43 namespace help
  44 {
  45
  46 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *//
  47
  48 template<class labelListType, class neiOp, class filterOp>
  49 void frontalMarking
  50 (
  51     labelListType& result,
  52     const label startingIndex,
  53     const neiOp& neighbourCalculator,
  54     const filterOp& selector
  55 )
  56 {
  57     //- add the starting element
  58     result.clear();
  59     result.append(startingIndex);
  60
  61     //- add the starting element to the front
  62     labelLongList front;
  63     front.append(startingIndex);
  64
  65     //- store information which element were already visited
  66     boolList alreadySelected(neighbourCalculator.size(), false);
  67
  68     //- start with frontal marking
  69     while( front.size() )
  70     {
  71         const label eLabel = front.removeLastElement();
  72
  73         //- find neighbours of the current element
  74         DynList<label> neighbours;
  75         neighbourCalculator(eLabel, neighbours);
  76
  77         forAll(neighbours, neiI)
  78         {
  79             const label nei = neighbours[neiI];
  80
  81             if( nei < 0 )
  82                 continue;
  83             if( alreadySelected[nei] )
  84                 continue;
  85
  86             if( selector(nei) )
  87             {
  88                 alreadySelected[nei] = true;
  89                 front.append(nei);
  90                 result.append(nei);
  91             }
  92         }
  93     }
  94 }
  95
  96 class graphNeiOp
  97 {
  98     // Private data
  99         //- const reference to VRWGraph
 100         const VRWGraph& neiGroups_;
 101
 102 public:
 103
 104     // Constructors
 105         //- Construct from VRWGraph
 106         inline graphNeiOp(const VRWGraph& neiGroups)
 107         :
 108             neiGroups_(neiGroups)
 109         {}
 110
 111     // Public member functions
 112         //- return the size of the graph
 113         inline label size() const
 114         {
 115             return neiGroups_.size();
 116         }
 117
 118         //- operator for finding neighbours
 119         inline void operator()(const label groupI, DynList<label>& ng) const
 120         {
 121             ng.setSize(neiGroups_.sizeOfRow(groupI));
 122             forAllRow(neiGroups_, groupI, i)
 123                 ng[i] = neiGroups_(groupI, i);
 124         }
 125 };
 126
 127 class graphSelectorOp
 128 {
 129     // Private data
 130         //- const reference to VRWGraph
 131         const VRWGraph& neiGroups_;
 132
 133 public:
 134
 135     // Constructors
 136         //- Construct from VRWGraph
 137         inline graphSelectorOp(const VRWGraph& neiGroups)
 138         :
 139             neiGroups_(neiGroups)
 140         {}
 141
 142     // Public member functions
 143         //- operator for selecting elements
 144         inline bool operator()(const label groupI) const
 145         {
 146             if( (groupI < 0) || (groupI >= neiGroups_.size()) )
 147                 return false;
 148
 149             return true;
 150         }
 151 };
 152
 153 template<class labelListType, class neiOp, class filterOp>
 154 label groupMarking
 155 (
 156     labelListType& elementInGroup,
 157     const neiOp& neighbourCalculator,
 158     const filterOp& selector
 159 )
 160 {
 161     label nGroups(0);
 162
 163     elementInGroup.setSize(neighbourCalculator.size());
 164     elementInGroup = -1;
 165
 166     VRWGraph neighbouringGroups;
 167
 168     label nThreads(1);
 169
 170     # ifdef USE_OMP
 171     //nThreads = 3 * omp_get_num_procs();
 172     # endif
 173
 174     DynList<label> nGroupsAtThread(nThreads, 0);
 175
 176     # ifdef USE_OMP
 177     # pragma omp parallel num_threads(nThreads)
 178     # endif
 179     {
 180         const label chunkSize =
 181             Foam::max(1, neighbourCalculator.size() / nThreads);
 182
 183         # ifdef USE_OMP
 184         const label threadI = omp_get_thread_num();
 185         # else
 186         const label threadI(0);
 187         # endif
 188
 189         LongList<std::pair<label, label> > threadCommPairs;
 190
 191         const label minEl = threadI * chunkSize;
 192
 193         label maxEl = minEl + chunkSize;
 194         if( threadI == (nThreads - 1) )
 195             maxEl = Foam::max(maxEl, neighbourCalculator.size());
 196
 197         label& groupI = nGroupsAtThread[threadI];
 198         groupI = 0;
 199
 200         for(label elI=minEl;elI<maxEl;++elI)
 201         {
 202             if( elementInGroup[elI] != -1 )
 203                 continue;
 204             if( !selector(elI) )
 205                 continue;
 206
 207             elementInGroup[elI] = groupI;
 208             labelLongList front;
 209             front.append(elI);
 210
 211             while( front.size() )
 212             {
 213                 const label eLabel = front.removeLastElement();
 214
 215                 DynList<label> neighbours;
 216                 neighbourCalculator(eLabel, neighbours);
 217
 218                 forAll(neighbours, neiI)
 219                 {
 220                     const label nei = neighbours[neiI];
 221
 222                     if( (nei < 0) || (elementInGroup[nei] != -1) )
 223                         continue;
 224
 225                     if( (nei < minEl) || (nei >= maxEl) )
 226                     {
 227                         //- this is a communication interface between
 228                         //- two threads
 229                         threadCommPairs.append(std::make_pair(elI, nei));
 230                     }
 231                     else if( selector(nei) )
 232                     {
 233                         //- this element is part of the same group
 234                         elementInGroup[nei] = groupI;
 235                         front.append(nei);
 236                     }
 237                 }
 238             }
 239
 240             ++groupI;
 241         }
 242
 243         # ifdef USE_OMP
 244         # pragma omp barrier
 245
 246         # pragma omp master
 247         {
 248             forAll(nGroupsAtThread, i)
 249                 nGroups += nGroupsAtThread[i];
 250         }
 251         # else
 252         nGroups = groupI;
 253         # endif
 254
 255         label startGroup(0);
 256         for(label i=0;i<threadI;++i)
 257             startGroup += nGroupsAtThread[i];
 258
 259         for(label elI=minEl;elI<maxEl;++elI)
 260             elementInGroup[elI] += startGroup;
 261
 262         # ifdef USE_OMP
 263         # pragma omp barrier
 264         # endif
 265
 266         //- find group to neighbouring groups addressing
 267         List<DynList<label> > localNeiGroups(nGroups);
 268         forAll(threadCommPairs, cfI)
 269         {
 270             const std::pair<label, label>& lp = threadCommPairs[cfI];
 271             const label groupI = elementInGroup[lp.first];
 272             const label neiGroup = elementInGroup[lp.second];
 273
 274             if( (neiGroup >= nGroups) || (groupI >= nGroups) )
 275                 FatalError << "neiGroup " << neiGroup
 276                     << " groupI " << groupI << " are >= than "
 277                     << "nGroups " << nGroups << abort(FatalError);
 278
 279             if( neiGroup != -1 )
 280             {
 281                 localNeiGroups[groupI].appendIfNotIn(neiGroup);
 282                 localNeiGroups[neiGroup].appendIfNotIn(groupI);
 283             }
 284         }
 285
 286         # ifdef USE_OMP
 287         # pragma omp critical
 288         # endif
 289         {
 290             neighbouringGroups.setSize(nGroups);
 291
 292             forAll(localNeiGroups, groupI)
 293             {
 294                 const DynList<label>& lGroups = localNeiGroups[groupI];
 295
 296                 neighbouringGroups.appendIfNotIn(groupI, groupI);
 297
 298                 forAll(lGroups, i)
 299                     neighbouringGroups.appendIfNotIn(groupI, lGroups[i]);
 300             }
 301         }
 302     }
 303
 304     forAll(neighbouringGroups, i)
 305     {
 306         labelList helper(neighbouringGroups.sizeOfRow(i));
 307         forAllRow(neighbouringGroups, i, j)
 308             helper[j] = neighbouringGroups(i, j);
 309
 310         sort(helper);
 311
 312         neighbouringGroups[i] = helper;
 313     }
 314
 315     //- start processing connections between the group and merge the connected
 316     //- ones into a new group
 317     DynList<label> globalGroupLabel;
 318     globalGroupLabel.setSize(nGroups);
 319     globalGroupLabel = -1;
 320
 321     //- reduce the information about the groups
 322     label counter(0);
 323
 324     forAll(neighbouringGroups, groupI)
 325     {
 326         if( globalGroupLabel[groupI] != -1 )
 327             continue;
 328
 329         DynList<label> connectedGroups;
 330         frontalMarking
 331         (
 332             connectedGroups,
 333             groupI,
 334             graphNeiOp(neighbouringGroups),
 335             graphSelectorOp(neighbouringGroups)
 336         );
 337
 338         forAll(connectedGroups, gI)
 339             globalGroupLabel[connectedGroups[gI]] = counter;
 340
 341         ++counter;
 342     }
 343
 344     nGroups = counter;
 345
 346     forAll(neighbouringGroups, groupI)
 347     {
 348         if( globalGroupLabel[groupI] != -1 )
 349             continue;
 350
 351         forAllRow(neighbouringGroups, groupI, ngI)
 352             globalGroupLabel[neighbouringGroups(groupI, ngI)] = counter;
 353
 354         ++counter;
 355     }
 356
 357     if( Pstream::parRun() )
 358     {
 359         //- reduce the groups over processors of an MPI run
 360         //- count the total number of groups over all processors
 361         labelList nGroupsAtProc(Pstream::nProcs());
 362         nGroupsAtProc[Pstream::myProcNo()] = nGroups;
 363
 364         Pstream::gatherList(nGroupsAtProc);
 365         Pstream::scatterList(nGroupsAtProc);
 366
 367         label startGroup(0), totalNumGroups(0);
 368         for(label procI=0;procI<Pstream::nProcs();++procI)
 369         {
 370             totalNumGroups += nGroupsAtProc[procI];
 371
 372             if( procI < Pstream::myProcNo() )
 373                 startGroup += nGroupsAtProc[procI];
 374         }
 375
 376         //- translate group labels
 377         forAll(globalGroupLabel, groupI)
 378             globalGroupLabel[groupI] += startGroup;
 379
 380         //- find the neighbouring groups
 381         //- collect groups on other processors
 382         //- this operator implements the algorithm for exchanging data
 383         //- over processors and collects information which groups
 384         //- are connected over inter-processor boundaries
 385         std::map<label, DynList<label> > neiGroups;
 386
 387         neighbourCalculator.collectGroups
 388         (
 389             neiGroups,
 390             elementInGroup,
 391             globalGroupLabel
 392         );
 393
 394         //- create a graph of connections
 395         List<List<labelPair> > globalNeiGroups(Pstream::nProcs());
 396
 397         DynList<labelPair> connsAtProc;
 398         for
 399         (
 400             std::map<label, DynList<label> >::const_iterator it =
 401             neiGroups.begin();
 402             it!=neiGroups.end();
 403             ++it
 404         )
 405         {
 406             const DynList<label>& ng = it->second;
 407
 408             forAll(ng, i)
 409                 connsAtProc.append(labelPair(it->first, ng[i]));
 410         }
 411
 412         //- copy the connections into the global neighbour list
 413         globalNeiGroups[Pstream::myProcNo()].setSize(connsAtProc.size());
 414
 415         forAll(connsAtProc, i)
 416             globalNeiGroups[Pstream::myProcNo()][i] = connsAtProc[i];
 417
 418         //- communicate partial graphs to the master processor
 419         Pstream::gatherList(globalNeiGroups);
 420
 421         labelList allGroupsNewLabel;
 422         if( Pstream::master() )
 423         {
 424             //- collect the graph of connections for the whole system
 425             VRWGraph allGroups(totalNumGroups);
 426             forAll(allGroups, i)
 427                 allGroups[i].append(i);
 428
 429             forAll(globalNeiGroups, procI)
 430             {
 431                 const List<labelPair>& connections = globalNeiGroups[procI];
 432
 433                 forAll(connections, i)
 434                 {
 435                     const labelPair& lp = connections[i];
 436
 437                     allGroups.appendIfNotIn(lp.first(), lp.second());
 438                     allGroups.appendIfNotIn(lp.second(), lp.first());
 439                 }
 440             }
 441
 442             //- assign a global label to each group
 443             allGroupsNewLabel.setSize(totalNumGroups);
 444             allGroupsNewLabel = -1;
 445             counter = 0;
 446
 447             forAll(allGroups, groupI)
 448             {
 449                 if( allGroupsNewLabel[groupI] != -1 )
 450                     continue;
 451
 452                 DynList<label> connectedGroups;
 453                 frontalMarking
 454                 (
 455                     connectedGroups,
 456                     groupI,
 457                     graphNeiOp(allGroups),
 458                     graphSelectorOp(allGroups)
 459                 );
 460
 461                 forAll(connectedGroups, gI)
 462                     allGroupsNewLabel[connectedGroups[gI]] = counter;
 463
 464                 ++counter;
 465             }
 466
 467             nGroups = counter;
 468         }
 469
 470         //- broadcast group labels from the master to other processors
 471         Pstream::scatter(nGroups);
 472         Pstream::scatter(allGroupsNewLabel);
 473
 474         //- assign correct group labels
 475         forAll(globalGroupLabel, groupI)
 476         {
 477             const label ngI = globalGroupLabel[groupI];
 478             globalGroupLabel[groupI] = allGroupsNewLabel[ngI];
 479         }
 480     }
 481
 482     //- set the global group label
 483     # ifdef USE_OMP
 484     # pragma omp parallel for schedule(dynamic, 50)
 485     # endif
 486     forAll(elementInGroup, elI)
 487     {
 488         if( elementInGroup[elI] < 0 )
 489             continue;
 490
 491         elementInGroup[elI] = globalGroupLabel[elementInGroup[elI]];
 492     }
 493
 494     return nGroups;
 495 }
 496
 497 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *//
 498
 499 } // End namespace help
 500
 501 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *//
 502
 503 } // End namespace Foam
 504
 505 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //