src/mesh/cfMesh/utilities/meshes/polyMeshGen/polyMeshGenCells.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 "polyMeshGenCells.H"
  29 #include "polyMeshGenAddressing.H"
  30 #include "IOobjectList.H"
  31 #include "cellSet.H"
  32 #include "demandDrivenData.H"
  33
  34 #include "labelPair.H"
  35
  36 # ifdef USE_OMP
  37 #include <omp.h>
  38 # endif
  39
  40 namespace Foam
  41 {
  42
  43 // * * * * * * * * * * Private member functions * * * * * * * * * * * * * * * //
  44
  45 void polyMeshGenCells::calculateOwnersAndNeighbours() const
  46 {
  47     if( ownerPtr_ || neighbourPtr_ )
  48         FatalErrorIn
  49         (
  50             "void polyMeshGenCells::calculateOwnersAndNeighbours() const"
  51         ) << "Owners and neighbours are already allocated" << abort(FatalError);
  52
  53     //- allocate owners
  54     ownerPtr_ =
  55         new labelIOList
  56         (
  57             IOobject
  58             (
  59                 "owner",
  60                 runTime_.constant(),
  61                 "polyMesh",
  62                 runTime_
  63             ),
  64             faces_.size()
  65         );
  66     labelIOList& own = *ownerPtr_;
  67
  68     //- allocate neighbours
  69     neighbourPtr_ =
  70         new labelIOList
  71         (
  72             IOobject
  73             (
  74                 "neighbour",
  75                 runTime_.constant(),
  76                 "polyMesh",
  77                 runTime_
  78             ),
  79             faces_.size()
  80         );
  81     labelIOList& nei = *neighbourPtr_;
  82
  83     //- start calculating owners and neighbours
  84     nIntFaces_ = 0;
  85
  86     # ifdef USE_OMP
  87     const label nThreads = 3 * omp_get_num_procs();
  88     const label chunkSize = faces_.size() / nThreads + 1;
  89     # else
  90     const label nThreads = 1;
  91     const label chunkSize = faces_.size();
  92     # endif
  93
  94     label nInternalFaces(0);
  95
  96     List<List<LongList<labelPair> > > dataForOtherThreads(nThreads);
  97
  98     # ifdef USE_OMP
  99     # pragma omp parallel num_threads(nThreads) reduction(+ : nInternalFaces)
 100     # endif
 101     {
 102         # ifdef USE_OMP
 103         const label threadI = omp_get_thread_num();
 104         # else
 105         const label threadI(0);
 106         # endif
 107
 108         const label startingFace = threadI * chunkSize;
 109         const label endFace =
 110             Foam::min(startingFace + chunkSize, faces_.size());
 111
 112         List<LongList<labelPair> >& dot = dataForOtherThreads[threadI];
 113         dot.setSize(nThreads);
 114
 115         for(label faceI=startingFace;faceI<endFace;++faceI)
 116         {
 117             own[faceI] = -1;
 118             nei[faceI] = -1;
 119         }
 120
 121         # ifdef USE_OMP
 122         # pragma omp for schedule(static)
 123         # endif
 124         forAll(cells_, cellI)
 125         {
 126             const cell& c = cells_[cellI];
 127
 128             forAll(c, fI)
 129             {
 130                 const label faceI = c[fI];
 131
 132                 const label threadNo = faceI / chunkSize;
 133
 134                 if( threadNo == threadI )
 135                 {
 136                     if( own[faceI] == -1 )
 137                     {
 138                         own[faceI] = cellI;
 139                     }
 140                     else if( nei[faceI] == -1 )
 141                     {
 142                         nei[faceI] = cellI;
 143                         ++nInternalFaces;
 144                     }
 145                     else
 146                     {
 147                         Serr << "Face " << faces_[faceI] << endl;
 148                         Serr << "Owner " << own[faceI] << endl;
 149                         Serr << "Neighbour " << nei[faceI] << endl;
 150                         Serr << "Current cell " << cellI << endl;
 151                         FatalErrorIn
 152                         (
 153                             "void polyMeshGenCells::"
 154                             "calculateOwnersAndNeighbours()"
 155                         ) << Pstream::myProcNo() << "Face " << faceI
 156                             << " appears in more than 2 cells!!"
 157                             << abort(FatalError);
 158                     }
 159                 }
 160                 else
 161                 {
 162                     dot[threadNo].append(labelPair(faceI, cellI));
 163                 }
 164             }
 165         }
 166
 167         # ifdef USE_OMP
 168         # pragma omp barrier
 169
 170         # pragma omp critical
 171         # endif
 172         for(label i=0;i<nThreads;++i)
 173         {
 174             const LongList<labelPair>& data =
 175                 dataForOtherThreads[i][threadI];
 176
 177             forAll(data, j)
 178             {
 179                 const label faceI = data[j].first();
 180                 const label cellI = data[j].second();
 181
 182                 if( own[faceI] == -1 )
 183                 {
 184                     own[faceI] = cellI;
 185                 }
 186                 else if( own[faceI] > cellI )
 187                 {
 188                     if( nei[faceI] == -1 )
 189                     {
 190                         nei[faceI] = own[faceI];
 191                         own[faceI] = cellI;
 192                         ++nInternalFaces;
 193                     }
 194                     else
 195                     {
 196                         Serr << "Face " << faces_[faceI] << endl;
 197                         Serr << "Owner " << own[faceI] << endl;
 198                         Serr << "Neighbour " << nei[faceI] << endl;
 199                         Serr << "Current cell " << cellI << endl;
 200                         FatalErrorIn
 201                         (
 202                             "void polyMeshGenCells::"
 203                             "calculateOwnersAndNeighbours()"
 204                         ) << Pstream::myProcNo() << "Face " << faceI
 205                             << " appears in more than 2 cells!!"
 206                             << abort(FatalError);
 207                     }
 208                 }
 209                 else if( nei[faceI] == -1 )
 210                 {
 211                     nei[faceI] = cellI;
 212                     ++nInternalFaces;
 213                 }
 214                 else
 215                 {
 216                     Serr << "Face " << faces_[faceI] << endl;
 217                     Serr << "Owner " << own[faceI] << endl;
 218                     Serr << "Neighbour " << nei[faceI] << endl;
 219                     Serr << "Current cell " << cellI << endl;
 220                     FatalErrorIn
 221                     (
 222                         "void polyMeshGenCells::"
 223                         "calculateOwnersAndNeighbours()"
 224                     ) << Pstream::myProcNo() << "Face " << faceI
 225                         << " appears in more than 2 cells!!"
 226                         << abort(FatalError);
 227                 }
 228             }
 229         }
 230     }
 231
 232     nIntFaces_ = nInternalFaces;
 233 }
 234
 235 void polyMeshGenCells::calculateAddressingData() const
 236 {
 237     if( !ownerPtr_ || !neighbourPtr_ )
 238     {
 239         # ifdef USE_OMP
 240         if( omp_in_parallel() )
 241             FatalErrorIn
 242             (
 243                 "inline label polyMeshGenCells::calculateAddressingData() const"
 244             ) << "Calculating addressing inside a parallel region."
 245                 << " This is not thread safe" << exit(FatalError);
 246         # endif
 247
 248         calculateOwnersAndNeighbours();
 249     }
 250
 251     addressingDataPtr_ = new polyMeshGenAddressing(*this);
 252 }
 253
 254 void polyMeshGenCells::clearOut() const
 255 {
 256     polyMeshGenFaces::clearOut();
 257     deleteDemandDrivenData(addressingDataPtr_);
 258 }
 259
 260 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 261 // Constructors
 262 //- Null constructor
 263 polyMeshGenCells::polyMeshGenCells(const Time& runTime)
 264 :
 265     polyMeshGenFaces(runTime),
 266     cells_(),
 267     cellSubsets_(),
 268     addressingDataPtr_(NULL)
 269 {
 270 }
 271
 272 //- Construct from components without the boundary
 273 polyMeshGenCells::polyMeshGenCells
 274 (
 275     const Time& runTime,
 276     const pointField& points,
 277     const faceList& faces,
 278     const cellList& cells
 279 )
 280 :
 281     polyMeshGenFaces(runTime, points, faces),
 282     cells_(),
 283     cellSubsets_(),
 284     addressingDataPtr_(NULL)
 285 {
 286     cells_ = cells;
 287 }
 288
 289 //- Construct from components with the boundary
 290 polyMeshGenCells::polyMeshGenCells
 291 (
 292     const Time& runTime,
 293     const pointField& points,
 294     const faceList& faces,
 295     const cellList& cells,
 296     const wordList& patchNames,
 297     const labelList& patchStart,
 298     const labelList& nFacesInPatch
 299 )
 300 :
 301     polyMeshGenFaces
 302     (
 303         runTime,
 304         points,
 305         faces,
 306         patchNames,
 307         patchStart,
 308         nFacesInPatch
 309     ),
 310     cells_(),
 311     cellSubsets_(),
 312     addressingDataPtr_(NULL)
 313 {
 314     cells_ = cells;
 315 }
 316
 317 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 318 // Destructor
 319 polyMeshGenCells::~polyMeshGenCells()
 320 {
 321     clearOut();
 322 }
 323
 324 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 325
 326 //- return addressing which may be needed
 327 const polyMeshGenAddressing& polyMeshGenCells::addressingData() const
 328 {
 329     if( !addressingDataPtr_ )
 330     {
 331         # ifdef USE_OMP
 332         if( omp_in_parallel() )
 333             FatalErrorIn
 334             (
 335                 "inline label polyMeshGenCells::addressingData() const"
 336             ) << "Calculating addressing inside a parallel region."
 337                 << " This is not thread safe" << exit(FatalError);
 338         # endif
 339
 340         calculateAddressingData();
 341     }
 342
 343     return *addressingDataPtr_;
 344 }
 345
 346 void polyMeshGenCells::clearAddressingData() const
 347 {
 348     deleteDemandDrivenData(addressingDataPtr_);
 349 }
 350
 351 label polyMeshGenCells::addCellSubset(const word& selName)
 352 {
 353     label id = cellSubsetIndex(selName);
 354     if( id >= 0 )
 355     {
 356         Warning << "Cell subset " << selName << " already exists!" << endl;
 357         return id;
 358     }
 359
 360     id = 0;
 361     for
 362     (
 363         std::map<label, meshSubset>::const_iterator it=cellSubsets_.begin();
 364         it!=cellSubsets_.end();
 365         ++it
 366     )
 367         id = Foam::max(id, it->first+1);
 368
 369     cellSubsets_.insert
 370     (
 371         std::make_pair
 372         (
 373             id,
 374             meshSubset(selName, meshSubset::CELLSUBSET)
 375         )
 376     );
 377
 378     return id;
 379 }
 380
 381 void polyMeshGenCells::removeCellSubset(const label setI)
 382 {
 383     if( cellSubsets_.find(setI) == cellSubsets_.end() )
 384         return;
 385
 386     cellSubsets_.erase(setI);
 387 }
 388
 389 word polyMeshGenCells::cellSubsetName(const label setI) const
 390 {
 391     std::map<label, meshSubset>::const_iterator it =
 392         cellSubsets_.find(setI);
 393     if( it == cellSubsets_.end() )
 394     {
 395         Warning << "Subset " << setI << " is not a cell subset" << endl;
 396         return word();
 397     }
 398
 399     return it->second.name();
 400 }
 401
 402 label polyMeshGenCells::cellSubsetIndex(const word& selName) const
 403 {
 404     std::map<label, meshSubset>::const_iterator it;
 405     for(it=cellSubsets_.begin();it!=cellSubsets_.end();++it)
 406     {
 407         if( it->second.name() == selName )
 408             return it->first;
 409     }
 410
 411     return -1;
 412 }
 413
 414 void polyMeshGenCells::read()
 415 {
 416     polyMeshGenFaces::read();
 417
 418     Info << "Starting creating cells" << endl;
 419     //- count the number of cells and create the cells
 420     label nCells(0);
 421     const labelList& own = this->owner();
 422     const labelList& nei = this->neighbour();
 423
 424     forAll(own, faceI)
 425     {
 426         if( own[faceI] >= nCells )
 427             nCells = own[faceI] + 1;
 428
 429         if( nei[faceI] >= nCells )
 430             nCells = nei[faceI] + 1;
 431     }
 432
 433     List<direction> nFacesInCell(nCells, direction(0));
 434     forAll(own, faceI)
 435         ++nFacesInCell[own[faceI]];
 436
 437     forAll(nei, faceI)
 438         if( nei[faceI] != -1 )
 439             ++nFacesInCell[nei[faceI]];
 440
 441     cells_.setSize(nCells);
 442     forAll(cells_, cellI)
 443         cells_[cellI].setSize(nFacesInCell[cellI]);
 444
 445     nFacesInCell = 0;
 446     forAll(own, faceI)
 447     {
 448         cells_[own[faceI]][nFacesInCell[own[faceI]]++] = faceI;
 449         if( nei[faceI] != -1 )
 450             cells_[nei[faceI]][nFacesInCell[nei[faceI]]++] = faceI;
 451     }
 452
 453     // read cell subsets
 454     IOobjectList allSets
 455     (
 456         runTime_,
 457         runTime_.constant(),
 458         "polyMesh/sets"
 459     );
 460
 461     wordList setNames = allSets.names("cellSet");
 462     forAll(setNames, setI)
 463     {
 464         IOobject* obj = allSets.lookup(setNames[setI]);
 465
 466         cellSet cSet(*obj);
 467
 468         const labelList content = cSet.toc();
 469         const label id = addCellSubset(setNames[setI]);
 470
 471         cellSubsets_[id].updateSubset(content);
 472     }
 473 }
 474
 475 void polyMeshGenCells::write() const
 476 {
 477     polyMeshGenFaces::write();
 478
 479     //- write cell subsets
 480     std::map<label, meshSubset>::const_iterator setIt;
 481     for(setIt=cellSubsets_.begin();setIt!=cellSubsets_.end();++setIt)
 482     {
 483         cellSet set
 484         (
 485             IOobject
 486             (
 487                 setIt->second.name(),
 488                 runTime_.constant(),
 489                 "polyMesh/sets",
 490                 runTime_,
 491                 IOobject::NO_READ,
 492                 IOobject::AUTO_WRITE
 493             )
 494         );
 495
 496         labelLongList containedElements;
 497         setIt->second.containedElements(containedElements);
 498
 499         forAll(containedElements, i)
 500             set.insert(containedElements[i]);
 501         set.write();
 502     }
 503 }
 504
 505 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 506
 507 } // End namespace Foam
 508
 509 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //