src/finiteVolume/fvMesh/singleCellFvMesh/singleCellFvMesh.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "singleCellFvMesh.H"
  27 #include "syncTools.H"
  28 #include "uindirectPrimitivePatch.H"
  29
  30 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  31
  32 // Conversion is a two step process:
  33 // - from original (fine) patch faces to agglomerations (aggloms might not
  34 //   be in correct patch order)
  35 // - from agglomerations to coarse patch faces
  36 void Foam::singleCellFvMesh::agglomerateMesh
  37 (
  38     const fvMesh& mesh,
  39     const labelListList& agglom
  40 )
  41 {
  42     const polyBoundaryMesh& oldPatches = mesh.boundaryMesh();
  43
  44     // Check agglomeration within patch face range and continuous
  45     labelList nAgglom(oldPatches.size(), 0);
  46
  47     forAll(oldPatches, patchI)
  48     {
  49         const polyPatch& pp = oldPatches[patchI];
  50         if (pp.size() > 0)
  51         {
  52             nAgglom[patchI] = max(agglom[patchI])+1;
  53
  54             forAll(pp, i)
  55             {
  56                 if (agglom[patchI][i] < 0  || agglom[patchI][i] >= pp.size())
  57                 {
  58                     FatalErrorIn
  59                     (
  60                         "singleCellFvMesh::agglomerateMesh(..)"
  61                     )   << "agglomeration on patch " << patchI
  62                         << " is out of range 0.." << pp.size()-1
  63                         << exit(FatalError);
  64                 }
  65             }
  66         }
  67     }
  68
  69     // Check agglomeration is sync
  70     {
  71         // Get neighbouring agglomeration
  72         labelList nbrAgglom(mesh.nFaces()-mesh.nInternalFaces());
  73         forAll(oldPatches, patchI)
  74         {
  75             const polyPatch& pp = oldPatches[patchI];
  76
  77             if (pp.coupled())
  78             {
  79                 label offset = pp.start()-mesh.nInternalFaces();
  80                 forAll(pp, i)
  81                 {
  82                     nbrAgglom[offset+i] = agglom[patchI][i];
  83                 }
  84             }
  85         }
  86         syncTools::swapBoundaryFaceList(mesh, nbrAgglom);
  87
  88
  89         // Get correspondence between this agglomeration and remote one
  90         Map<label> localToNbr(nbrAgglom.size()/10);
  91
  92         forAll(oldPatches, patchI)
  93         {
  94             const polyPatch& pp = oldPatches[patchI];
  95
  96             if (pp.coupled())
  97             {
  98                 label offset = pp.start()-mesh.nInternalFaces();
  99
 100                 forAll(pp, i)
 101                 {
 102                     label bFaceI = offset+i;
 103                     label myZone = agglom[patchI][i];
 104                     label nbrZone = nbrAgglom[bFaceI];
 105
 106                     Map<label>::const_iterator iter = localToNbr.find(myZone);
 107
 108                     if (iter == localToNbr.end())
 109                     {
 110                         // First occurence of this zone. Store correspondence
 111                         // to remote zone number.
 112                         localToNbr.insert(myZone, nbrZone);
 113                     }
 114                     else
 115                     {
 116                         // Check that zone numbers are still the same.
 117                         if (iter() != nbrZone)
 118                         {
 119                             FatalErrorIn
 120                             (
 121                                 "singleCellFvMesh::agglomerateMesh(..)"
 122                             )   << "agglomeration is not synchronised across"
 123                                 << " coupled patch " << pp.name()
 124                                 << endl
 125                                 << "Local agglomeration " << myZone
 126                                 << ". Remote agglomeration " << nbrZone
 127                                 << exit(FatalError);
 128                         }
 129                     }
 130                 }
 131             }
 132         }
 133     }
 134
 135
 136     label coarseI = 0;
 137     forAll(nAgglom, patchI)
 138     {
 139         coarseI += nAgglom[patchI];
 140     }
 141     // New faces
 142     faceList patchFaces(coarseI);
 143     // New patch start and size
 144     labelList patchStarts(oldPatches.size());
 145     labelList patchSizes(oldPatches.size());
 146
 147     // From new patch face back to agglomeration
 148     patchFaceMap_.setSize(oldPatches.size());
 149
 150     // From fine face to coarse face (or -1)
 151     reverseFaceMap_.setSize(mesh.nFaces());
 152     reverseFaceMap_.labelList::operator=(-1);
 153
 154     // Face counter
 155     coarseI = 0;
 156
 157
 158     forAll(oldPatches, patchI)
 159     {
 160         patchStarts[patchI] = coarseI;
 161
 162         const polyPatch& pp = oldPatches[patchI];
 163
 164         if (pp.size() > 0)
 165         {
 166             patchFaceMap_[patchI].setSize(nAgglom[patchI]);
 167
 168             // Patchfaces per agglomeration
 169             labelListList agglomToPatch
 170             (
 171                 invertOneToMany(nAgglom[patchI], agglom[patchI])
 172             );
 173
 174             // From agglomeration to compact patch face
 175             labelList agglomToFace(nAgglom[patchI], -1);
 176
 177             forAll(pp, i)
 178             {
 179                 label myAgglom = agglom[patchI][i];
 180
 181                 if (agglomToFace[myAgglom] == -1)
 182                 {
 183                     // Agglomeration not yet done. We now have:
 184                     // - coarseI                  : current coarse mesh face
 185                     // - patchStarts[patchI]      : coarse mesh patch start
 186                     // - myAgglom                 : agglomeration
 187                     // -  agglomToPatch[myAgglom] : fine mesh faces for zone
 188                     label coarsePatchFaceI = coarseI - patchStarts[patchI];
 189                     patchFaceMap_[patchI][coarsePatchFaceI] = myAgglom;
 190                     agglomToFace[myAgglom] = coarsePatchFaceI;
 191
 192                     const labelList& fineFaces = agglomToPatch[myAgglom];
 193
 194                     // Create overall map from fine mesh faces to coarseI.
 195                     forAll(fineFaces, fineI)
 196                     {
 197                         reverseFaceMap_[pp.start()+fineFaces[fineI]] = coarseI;
 198                     }
 199
 200                     // Construct single face
 201                     uindirectPrimitivePatch upp
 202                     (
 203                         UIndirectList<face>(pp, fineFaces),
 204                         pp.points()
 205                     );
 206
 207                     if (upp.edgeLoops().size() != 1)
 208                     {
 209                         FatalErrorIn
 210                         (
 211                             "singleCellFvMesh::agglomerateMesh(..)"
 212                         )   << "agglomeration does not create a"
 213                             << " single, non-manifold"
 214                             << " face for agglomeration " << myAgglom
 215                             << " on patch " <<  patchI
 216                             << exit(FatalError);
 217                     }
 218
 219                     patchFaces[coarseI++] = face
 220                     (
 221                         renumber
 222                         (
 223                             upp.meshPoints(),
 224                             upp.edgeLoops()[0]
 225                         )
 226                     );
 227                 }
 228             }
 229         }
 230
 231         patchSizes[patchI] = coarseI-patchStarts[patchI];
 232     }
 233
 234     //Pout<< "patchStarts:" << patchStarts << endl;
 235     //Pout<< "patchSizes:" << patchSizes << endl;
 236
 237     // Compact numbering for points
 238     reversePointMap_.setSize(mesh.nPoints());
 239     reversePointMap_.labelList::operator=(-1);
 240     label newI = 0;
 241
 242     forAll(patchFaces, coarseI)
 243     {
 244         face& f = patchFaces[coarseI];
 245
 246         forAll(f, fp)
 247         {
 248             if (reversePointMap_[f[fp]] == -1)
 249             {
 250                 reversePointMap_[f[fp]] = newI++;
 251             }
 252
 253             f[fp] = reversePointMap_[f[fp]];
 254         }
 255     }
 256
 257     pointMap_ = invert(newI, reversePointMap_);
 258
 259     // Subset used points
 260     pointField boundaryPoints(mesh.points(), pointMap_);
 261
 262     // Add patches (on still zero sized mesh)
 263     List<polyPatch*> newPatches(oldPatches.size());
 264     forAll(oldPatches, patchI)
 265     {
 266         newPatches[patchI] = oldPatches[patchI].clone
 267         (
 268             boundaryMesh(),
 269             patchI,
 270             0,
 271             0
 272         ).ptr();
 273     }
 274     addFvPatches(newPatches);
 275
 276     // Owner, neighbour is trivial
 277     labelList owner(patchFaces.size(), 0);
 278     labelList neighbour(0);
 279
 280
 281     // actually change the mesh
 282     resetPrimitives
 283     (
 284         xferMove(boundaryPoints),
 285         xferMove(patchFaces),
 286         xferMove(owner),
 287         xferMove(neighbour),
 288         patchSizes,
 289         patchStarts,
 290         true                //syncPar
 291     );
 292
 293
 294     // Adapt the zones
 295     cellZones().clear();
 296     cellZones().setSize(mesh.cellZones().size());
 297     {
 298         forAll(mesh.cellZones(), zoneI)
 299         {
 300             const cellZone& oldFz = mesh.cellZones()[zoneI];
 301
 302             DynamicList<label> newAddressing;
 303
 304             //Note: uncomment if you think it makes sense. Note that value
 305             // of cell0 is the average.
 306             //// Was old cell0 in this cellZone?
 307             //if (oldFz.localID(0) != -1)
 308             //{
 309             //    newAddressing.append(0);
 310             //}
 311
 312             cellZones().set
 313             (
 314                 zoneI,
 315                 oldFz.clone
 316                 (
 317                     newAddressing,
 318                     zoneI,
 319                     cellZones()
 320                 )
 321             );
 322         }
 323     }
 324
 325     faceZones().clear();
 326     faceZones().setSize(mesh.faceZones().size());
 327     {
 328         forAll(mesh.faceZones(), zoneI)
 329         {
 330             const faceZone& oldFz = mesh.faceZones()[zoneI];
 331
 332             DynamicList<label> newAddressing(oldFz.size());
 333             DynamicList<bool> newFlipMap(oldFz.size());
 334
 335             forAll(oldFz, i)
 336             {
 337                 label newFaceI = reverseFaceMap_[oldFz[i]];
 338
 339                 if (newFaceI != -1)
 340                 {
 341                     newAddressing.append(newFaceI);
 342                     newFlipMap.append(oldFz.flipMap()[i]);
 343                 }
 344             }
 345
 346             faceZones().set
 347             (
 348                 zoneI,
 349                 oldFz.clone
 350                 (
 351                     newAddressing,
 352                     newFlipMap,
 353                     zoneI,
 354                     faceZones()
 355                 )
 356             );
 357         }
 358     }
 359
 360
 361     pointZones().clear();
 362     pointZones().setSize(mesh.pointZones().size());
 363     {
 364         forAll(mesh.pointZones(), zoneI)
 365         {
 366             const pointZone& oldFz = mesh.pointZones()[zoneI];
 367
 368             DynamicList<label> newAddressing(oldFz.size());
 369
 370             forAll(oldFz, i)
 371             {
 372                 label newPointI  = reversePointMap_[oldFz[i]];
 373                 if (newPointI != -1)
 374                 {
 375                     newAddressing.append(newPointI);
 376                 }
 377             }
 378
 379             pointZones().set
 380             (
 381                 zoneI,
 382                 oldFz.clone
 383                 (
 384                     pointZones(),
 385                     zoneI,
 386                     newAddressing
 387                 )
 388             );
 389         }
 390     }
 391 }
 392
 393
 394 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 395
 396 Foam::singleCellFvMesh::singleCellFvMesh
 397 (
 398     const IOobject& io,
 399     const fvMesh& mesh
 400 )
 401 :
 402     fvMesh
 403     (
 404         io,
 405         xferCopy(pointField()), //points
 406         xferCopy(faceList()),   //faces
 407         xferCopy(labelList()),  //allOwner
 408         xferCopy(labelList()),  //allNeighbour
 409         false                   //syncPar
 410     ),
 411     patchFaceAgglomeration_
 412     (
 413         IOobject
 414         (
 415             "patchFaceAgglomeration",
 416             io.instance(),
 417             fvMesh::meshSubDir,
 418             *this,
 419             io.readOpt(),
 420             io.writeOpt()
 421         ),
 422         0
 423     ),
 424     patchFaceMap_
 425     (
 426         IOobject
 427         (
 428             "patchFaceMap",
 429             io.instance(),
 430             fvMesh::meshSubDir,
 431             *this,
 432             io.readOpt(),
 433             io.writeOpt()
 434         ),
 435         mesh.boundaryMesh().size()
 436     ),
 437     reverseFaceMap_
 438     (
 439         IOobject
 440         (
 441             "reverseFaceMap",
 442             io.instance(),
 443             fvMesh::meshSubDir,
 444             *this,
 445             io.readOpt(),
 446             io.writeOpt()
 447         ),
 448         mesh.nFaces()
 449     ),
 450     pointMap_
 451     (
 452         IOobject
 453         (
 454             "pointMap",
 455             io.instance(),
 456             fvMesh::meshSubDir,
 457             *this,
 458             io.readOpt(),
 459             io.writeOpt()
 460         ),
 461         mesh.nPoints()
 462     ),
 463     reversePointMap_
 464     (
 465         IOobject
 466         (
 467             "reversePointMap",
 468             io.instance(),
 469             fvMesh::meshSubDir,
 470             *this,
 471             io.readOpt(),
 472             io.writeOpt()
 473         ),
 474         mesh.nPoints()
 475     )
 476 {
 477     const polyBoundaryMesh& oldPatches = mesh.boundaryMesh();
 478
 479     labelListList agglom(oldPatches.size());
 480
 481     forAll(oldPatches, patchI)
 482     {
 483         agglom[patchI] = identity(oldPatches[patchI].size());
 484     }
 485
 486     agglomerateMesh(mesh, agglom);
 487 }
 488
 489
 490 Foam::singleCellFvMesh::singleCellFvMesh
 491 (
 492     const IOobject& io,
 493     const fvMesh& mesh,
 494     const labelListList& patchFaceAgglomeration
 495 )
 496 :
 497     fvMesh
 498     (
 499         io,
 500         xferCopy(pointField()), //points
 501         xferCopy(faceList()),   //faces
 502         xferCopy(labelList()),  //allOwner
 503         xferCopy(labelList()),  //allNeighbour
 504         false                   //syncPar
 505     ),
 506     patchFaceAgglomeration_
 507     (
 508         IOobject
 509         (
 510             "patchFaceAgglomeration",
 511             io.instance(),
 512             fvMesh::meshSubDir,
 513             *this,
 514             io.readOpt(),
 515             io.writeOpt()
 516         ),
 517         patchFaceAgglomeration
 518     ),
 519     patchFaceMap_
 520     (
 521         IOobject
 522         (
 523             "patchFaceMap",
 524             io.instance(),
 525             fvMesh::meshSubDir,
 526             *this,
 527             io.readOpt(),
 528             io.writeOpt()
 529         ),
 530         mesh.boundaryMesh().size()
 531     ),
 532     reverseFaceMap_
 533     (
 534         IOobject
 535         (
 536             "reverseFaceMap",
 537             io.instance(),
 538             fvMesh::meshSubDir,
 539             *this,
 540             io.readOpt(),
 541             io.writeOpt()
 542         ),
 543         mesh.nFaces()
 544     ),
 545     pointMap_
 546     (
 547         IOobject
 548         (
 549             "pointMap",
 550             io.instance(),
 551             fvMesh::meshSubDir,
 552             *this,
 553             io.readOpt(),
 554             io.writeOpt()
 555         ),
 556         mesh.nPoints()
 557     ),
 558     reversePointMap_
 559     (
 560         IOobject
 561         (
 562             "reversePointMap",
 563             io.instance(),
 564             fvMesh::meshSubDir,
 565             *this,
 566             io.readOpt(),
 567             io.writeOpt()
 568         ),
 569         mesh.nPoints()
 570     )
 571 {
 572     agglomerateMesh(mesh, patchFaceAgglomeration);
 573 }
 574
 575
 576 Foam::singleCellFvMesh::singleCellFvMesh(const IOobject& io)
 577 :
 578     fvMesh(io),
 579     patchFaceAgglomeration_
 580     (
 581         IOobject
 582         (
 583             "patchFaceAgglomeration",
 584             io.instance(),
 585             fvMesh::meshSubDir,
 586             *this,
 587             io.readOpt(),
 588             io.writeOpt()
 589         )
 590     ),
 591     patchFaceMap_
 592     (
 593         IOobject
 594         (
 595             "patchFaceMap",
 596             io.instance(),
 597             fvMesh::meshSubDir,
 598             *this,
 599             io.readOpt(),
 600             io.writeOpt()
 601         )
 602     ),
 603     reverseFaceMap_
 604     (
 605         IOobject
 606         (
 607             "reverseFaceMap",
 608             io.instance(),
 609             fvMesh::meshSubDir,
 610             *this,
 611             io.readOpt(),
 612             io.writeOpt()
 613         )
 614     ),
 615     pointMap_
 616     (
 617         IOobject
 618         (
 619             "pointMap",
 620             io.instance(),
 621             fvMesh::meshSubDir,
 622             *this,
 623             io.readOpt(),
 624             io.writeOpt()
 625         )
 626     ),
 627     reversePointMap_
 628     (
 629         IOobject
 630         (
 631             "reversePointMap",
 632             io.instance(),
 633             fvMesh::meshSubDir,
 634             *this,
 635             io.readOpt(),
 636             io.writeOpt()
 637         )
 638     )
 639 {}
 640
 641
 642 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 643
 644
 645
 646 // ************************************************************************* //