src/dynamicMesh/polyTopoChange/polyTopoChange/refinementHistory.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
   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 the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     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, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*---------------------------------------------------------------------------*/
  26
  27 #include "DynamicList.H"
  28 #include "refinementHistory.H"
  29 #include "ListOps.H"
  30 #include "mapPolyMesh.H"
  31 #include "mapDistributePolyMesh.H"
  32 #include "polyMesh.H"
  33
  34 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  35
  36 namespace Foam
  37 {
  38
  39 defineTypeNameAndDebug(refinementHistory, 0);
  40
  41 }
  42
  43
  44 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  45
  46 void Foam::refinementHistory::writeEntry
  47 (
  48     const List<splitCell8>& splitCells,
  49     const splitCell8& split
  50 )
  51 {
  52     // Write me:
  53     if (split.addedCellsPtr_.valid())
  54     {
  55         Pout<< "parent:" << split.parent_
  56             << " subCells:" << split.addedCellsPtr_()
  57             << endl;
  58     }
  59     else
  60     {
  61         Pout<< "parent:" << split.parent_
  62             << " no subcells"
  63             << endl;
  64     }
  65
  66     if (split.parent_ >= 0)
  67     {
  68         Pout<< "parent data:" << endl;
  69         // Write my parent
  70         string oldPrefix = Pout.prefix();
  71         Pout.prefix() = "  " + oldPrefix;
  72         writeEntry(splitCells, splitCells[split.parent_]);
  73         Pout.prefix() = oldPrefix;
  74     }
  75 }
  76
  77
  78 void Foam::refinementHistory::writeDebug
  79 (
  80     const labelList& visibleCells,
  81     const List<splitCell8>& splitCells
  82 )
  83 {
  84     string oldPrefix = Pout.prefix();
  85     Pout.prefix() = "";
  86
  87     forAll(visibleCells, cellI)
  88     {
  89         label index = visibleCells[cellI];
  90
  91         if (index >= 0)
  92         {
  93             Pout<< "Cell from refinement:" << cellI << " index:" << index
  94                 << endl;
  95
  96             string oldPrefix = Pout.prefix();
  97             Pout.prefix() = "  " + oldPrefix;
  98             writeEntry(splitCells, splitCells[index]);
  99             Pout.prefix() = oldPrefix;
 100         }
 101         else
 102         {
 103             Pout<< "Unrefined cell:" << cellI << " index:" << index << endl;
 104         }
 105     }
 106     Pout.prefix() = oldPrefix;
 107 }
 108
 109
 110 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 111
 112 //- Construct null
 113 Foam::refinementHistory::splitCell8::splitCell8()
 114 :
 115     parent_(-1),
 116     addedCellsPtr_(NULL)
 117 {}
 118
 119
 120 //- Construct as child element of parent
 121 Foam::refinementHistory::splitCell8::splitCell8(const label parent)
 122 :
 123     parent_(parent),
 124     addedCellsPtr_(NULL)
 125 {}
 126
 127
 128 //- Construct from Istream
 129 Foam::refinementHistory::splitCell8::splitCell8(Istream& is)
 130 {
 131     is >> *this;
 132 }
 133
 134
 135 //- Construct as (deep) copy.
 136 Foam::refinementHistory::splitCell8::splitCell8(const splitCell8& sc)
 137 :
 138     parent_(sc.parent_),
 139     addedCellsPtr_
 140     (
 141         sc.addedCellsPtr_.valid()
 142       ? new FixedList<label, 8>(sc.addedCellsPtr_())
 143       : NULL
 144     )
 145 {}
 146
 147
 148 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 149
 150 Foam::Istream& Foam::operator>>(Istream& is, refinementHistory::splitCell8& sc)
 151 {
 152     labelList addedCells;
 153
 154     is >> sc.parent_ >> addedCells;
 155
 156     if (addedCells.size())
 157     {
 158         sc.addedCellsPtr_.reset(new FixedList<label, 8>(addedCells));
 159     }
 160     else
 161     {
 162         sc.addedCellsPtr_.reset(NULL);
 163     }
 164
 165     return is;
 166 }
 167
 168
 169 Foam::Ostream& Foam::operator<<
 170 (
 171     Ostream& os,
 172     const refinementHistory::splitCell8& sc
 173 )
 174 {
 175     // Output as labelList so we can have 0 sized lists. Alternative is to
 176     // output as fixedlist with e.g. -1 elements and check for this upon
 177     // reading. However would cause much more data to be transferred.
 178
 179     if (sc.addedCellsPtr_.valid())
 180     {
 181         return os
 182             << sc.parent_
 183             << token::SPACE
 184             << labelList(sc.addedCellsPtr_());
 185     }
 186     else
 187     {
 188         return os << sc.parent_ << token::SPACE << labelList(0);
 189     }
 190 }
 191
 192
 193 void Foam::refinementHistory::checkIndices() const
 194 {
 195     // Check indices.
 196     forAll(visibleCells_, i)
 197     {
 198         if (visibleCells_[i] < 0 && visibleCells_[i] >= splitCells_.size())
 199         {
 200             FatalErrorIn("refinementHistory::checkIndices() const")
 201                 << "Illegal entry " << visibleCells_[i]
 202                 << " in visibleCells at location" << i << nl
 203                 << "It points outside the range of splitCells : 0.."
 204                 << splitCells_.size()-1
 205                 << abort(FatalError);
 206         }
 207     }
 208 }
 209
 210
 211 Foam::label Foam::refinementHistory::allocateSplitCell
 212 (
 213     const label parent,
 214     const label i
 215 )
 216 {
 217     label index = -1;
 218
 219     if (freeSplitCells_.size())
 220     {
 221         index = freeSplitCells_.remove();
 222
 223         splitCells_[index] = splitCell8(parent);
 224     }
 225     else
 226     {
 227         index = splitCells_.size();
 228
 229         splitCells_.append(splitCell8(parent));
 230     }
 231
 232
 233     // Update the parent field
 234     if (parent >= 0)
 235     {
 236         splitCell8& parentSplit = splitCells_[parent];
 237
 238         if (parentSplit.addedCellsPtr_.empty())
 239         {
 240             // Allocate storage on parent for the 8 subcells.
 241             parentSplit.addedCellsPtr_.reset(new FixedList<label, 8>(-1));
 242         }
 243
 244
 245         // Store me on my parent
 246         FixedList<label, 8>& parentSplits = parentSplit.addedCellsPtr_();
 247
 248         parentSplits[i] = index;
 249     }
 250
 251     return index;
 252 }
 253
 254
 255 void Foam::refinementHistory::freeSplitCell(const label index)
 256 {
 257     splitCell8& split = splitCells_[index];
 258
 259     // Make sure parent does not point to me anymore.
 260     if (split.parent_ >= 0)
 261     {
 262         autoPtr<FixedList<label, 8> >& subCellsPtr =
 263             splitCells_[split.parent_].addedCellsPtr_;
 264
 265         if (subCellsPtr.valid())
 266         {
 267             FixedList<label, 8>& subCells = subCellsPtr();
 268
 269             label myPos = findIndex(subCells, index);
 270
 271             if (myPos == -1)
 272             {
 273                 FatalErrorIn("refinementHistory::freeSplitCell")
 274                     << "Problem: cannot find myself in"
 275                     << " parents' children" << abort(FatalError);
 276             }
 277             else
 278             {
 279                 subCells[myPos] = -1;
 280             }
 281         }
 282     }
 283
 284     // Mark splitCell as free
 285     split.parent_ = -2;
 286
 287     // Add to cache of free splitCells
 288     freeSplitCells_.append(index);
 289 }
 290
 291
 292 // Mark entry in splitCells. Recursively mark its parent and subs.
 293 void Foam::refinementHistory::markSplit
 294 (
 295     const label index,
 296     labelList& oldToNew,
 297     DynamicList<splitCell8>& newSplitCells
 298 ) const
 299 {
 300     if (oldToNew[index] == -1)
 301     {
 302         // Not yet compacted.
 303
 304         const splitCell8& split = splitCells_[index];
 305
 306         oldToNew[index] = newSplitCells.size();
 307         newSplitCells.append(split);
 308
 309         if (split.parent_ >= 0)
 310         {
 311             markSplit(split.parent_, oldToNew, newSplitCells);
 312         }
 313         if (split.addedCellsPtr_.valid())
 314         {
 315             const FixedList<label, 8>& splits = split.addedCellsPtr_();
 316
 317             forAll(splits, i)
 318             {
 319                 if (splits[i] >= 0)
 320                 {
 321                     markSplit(splits[i], oldToNew, newSplitCells);
 322                 }
 323             }
 324         }
 325     }
 326 }
 327
 328
 329 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 330
 331 Foam::refinementHistory::refinementHistory(const IOobject& io)
 332 :
 333     regIOobject(io)
 334 {
 335     if
 336     (
 337         io.readOpt() == IOobject::MUST_READ
 338      || (io.readOpt() == IOobject::READ_IF_PRESENT && headerOk())
 339     )
 340     {
 341         readStream(typeName) >> *this;
 342         close();
 343     }
 344
 345     if (debug)
 346     {
 347         Pout<< "refinementHistory::refinementHistory :"
 348             << " constructed history from IOobject :"
 349             << " splitCells:" << splitCells_.size()
 350             << " visibleCells:" << visibleCells_.size()
 351             << endl;
 352     }
 353 }
 354
 355
 356 //- Read or construct
 357 Foam::refinementHistory::refinementHistory
 358 (
 359     const IOobject& io,
 360     const List<splitCell8>& splitCells,
 361     const labelList& visibleCells
 362 )
 363 :
 364     regIOobject(io),
 365     splitCells_(splitCells),
 366     freeSplitCells_(0),
 367     visibleCells_(visibleCells)
 368 {
 369     if
 370     (
 371         io.readOpt() == IOobject::MUST_READ
 372      || (io.readOpt() == IOobject::READ_IF_PRESENT && headerOk())
 373     )
 374     {
 375         readStream(typeName) >> *this;
 376         close();
 377     }
 378
 379     // Check indices.
 380     checkIndices();
 381
 382     if (debug)
 383     {
 384         Pout<< "refinementHistory::refinementHistory :"
 385             << " constructed history from IOobject or components :"
 386             << " splitCells:" << splitCells_.size()
 387             << " visibleCells:" << visibleCells_.size()
 388             << endl;
 389     }
 390 }
 391
 392
 393 // Construct from initial number of cells (all visible)
 394 Foam::refinementHistory::refinementHistory
 395 (
 396     const IOobject& io,
 397     const label nCells
 398 )
 399 :
 400     regIOobject(io),
 401     freeSplitCells_(0)
 402 {
 403     if
 404     (
 405         io.readOpt() == IOobject::MUST_READ
 406      || (io.readOpt() == IOobject::READ_IF_PRESENT && headerOk())
 407     )
 408     {
 409         readStream(typeName) >> *this;
 410         close();
 411     }
 412     else
 413     {
 414         visibleCells_.setSize(nCells);
 415         splitCells_.setCapacity(nCells);
 416
 417         for (label cellI = 0; cellI < nCells; cellI++)
 418         {
 419             visibleCells_[cellI] = cellI;
 420             splitCells_.append(splitCell8());
 421         }
 422     }
 423
 424     // Check indices.
 425     checkIndices();
 426
 427     if (debug)
 428     {
 429         Pout<< "refinementHistory::refinementHistory :"
 430             << " constructed history from IOobject or initial size :"
 431             << " splitCells:" << splitCells_.size()
 432             << " visibleCells:" << visibleCells_.size()
 433             << endl;
 434     }
 435 }
 436
 437
 438 // Construct as copy
 439 Foam::refinementHistory::refinementHistory
 440 (
 441     const IOobject& io,
 442     const refinementHistory& rh
 443 )
 444 :
 445     regIOobject(io),
 446     splitCells_(rh.splitCells()),
 447     freeSplitCells_(rh.freeSplitCells()),
 448     visibleCells_(rh.visibleCells())
 449 {
 450     if (debug)
 451     {
 452         Pout<< "refinementHistory::refinementHistory : constructed initial"
 453             << " history." << endl;
 454     }
 455 }
 456
 457
 458 // Construct from Istream
 459 Foam::refinementHistory::refinementHistory(const IOobject& io, Istream& is)
 460 :
 461     regIOobject(io),
 462     splitCells_(is),
 463     freeSplitCells_(0),
 464     visibleCells_(is)
 465 {
 466     // Check indices.
 467     checkIndices();
 468
 469     if (debug)
 470     {
 471         Pout<< "refinementHistory::refinementHistory :"
 472             << " constructed history from Istream"
 473             << " splitCells:" << splitCells_.size()
 474             << " visibleCells:" << visibleCells_.size()
 475             << endl;
 476     }
 477 }
 478
 479
 480 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 481
 482 void Foam::refinementHistory::resize(const label size)
 483 {
 484     label oldSize = visibleCells_.size();
 485
 486     if (debug)
 487     {
 488         Pout<< "refinementHistory::resize from " << oldSize << " to " << size
 489             << " cells" << endl;
 490     }
 491
 492     visibleCells_.setSize(size);
 493
 494     // Set additional elements to -1.
 495     for (label i = oldSize; i < visibleCells_.size(); i++)
 496     {
 497         visibleCells_[i] = -1;
 498     }
 499 }
 500
 501
 502 void Foam::refinementHistory::updateMesh(const mapPolyMesh& map)
 503 {
 504     if (active())
 505     {
 506         const labelList& reverseCellMap = map.reverseCellMap();
 507
 508         // Note that only the live cells need to be renumbered.
 509
 510         labelList newVisibleCells(map.cellMap().size(), -1);
 511
 512         forAll(visibleCells_, cellI)
 513         {
 514             if (visibleCells_[cellI] != -1)
 515             {
 516                 label index = visibleCells_[cellI];
 517
 518                 // Check not already set
 519                 if (splitCells_[index].addedCellsPtr_.valid())
 520                 {
 521                     FatalErrorIn
 522                     (
 523                         "refinementHistory::updateMesh(const mapPolyMesh&)"
 524                     )   << "Problem" << abort(FatalError);
 525                 }
 526
 527                 label newCellI = reverseCellMap[cellI];
 528
 529                 if (newCellI >= 0)
 530                 {
 531                     newVisibleCells[newCellI] = index;
 532                 }
 533             }
 534         }
 535
 536         if (debug)
 537         {
 538             Pout<< "refinementHistory::updateMesh : from "
 539                 << visibleCells_.size()
 540                 << " to " << newVisibleCells.size()
 541                 << " cells" << endl;
 542         }
 543
 544         visibleCells_.transfer(newVisibleCells);
 545     }
 546 }
 547
 548
 549 // Update numbering for subsetting
 550 void Foam::refinementHistory::subset
 551 (
 552     const labelList& pointMap,
 553     const labelList& faceMap,
 554     const labelList& cellMap
 555 )
 556 {
 557     if (active())
 558     {
 559         labelList newVisibleCells(cellMap.size(), -1);
 560
 561         forAll(newVisibleCells, cellI)
 562         {
 563             label oldCellI = cellMap[cellI];
 564
 565             label index = visibleCells_[oldCellI];
 566
 567             // Check that cell is live (so its parent has no refinement)
 568             if (index >= 0 && splitCells_[index].addedCellsPtr_.valid())
 569             {
 570                 FatalErrorIn
 571                 (
 572                     "refinementHistory::subset"
 573                     "(const labelList&, const labelList&, const labelList&)"
 574                 )   << "Problem" << abort(FatalError);
 575             }
 576
 577             newVisibleCells[cellI] = index;
 578         }
 579
 580         if (debug)
 581         {
 582             Pout<< "refinementHistory::updateMesh : from "
 583                 << visibleCells_.size()
 584                 << " to " << newVisibleCells.size()
 585                 << " cells" << endl;
 586         }
 587
 588         visibleCells_.transfer(newVisibleCells);
 589     }
 590 }
 591
 592
 593 void Foam::refinementHistory::countProc
 594 (
 595     const label index,
 596     const label newProcNo,
 597     labelList& splitCellProc,
 598     labelList& splitCellNum
 599 ) const
 600 {
 601     if (splitCellProc[index] != newProcNo)
 602     {
 603         // Different destination processor from other cells using this
 604         // parent. Reset count.
 605         splitCellProc[index] = newProcNo;
 606         splitCellNum[index] = 1;
 607     }
 608     else
 609     {
 610         splitCellNum[index]++;
 611
 612         // Increment parent if whole splitCell moves to same processor
 613         if (splitCellNum[index] == 8)
 614         {
 615             Pout<< "Moving " << splitCellNum[index]
 616                 << " cells originating from cell " << index
 617                 << " from processor " << Pstream::myProcNo()
 618                 << " to processor " << splitCellProc[index]
 619                 << endl;
 620
 621             label parent = splitCells_[index].parent_;
 622
 623             if (parent >= 0)
 624             {
 625                 string oldPrefix = Pout.prefix();
 626                 Pout.prefix() = "  " + oldPrefix;
 627
 628                 countProc(parent, newProcNo, splitCellProc, splitCellNum);
 629
 630                 Pout.prefix() = oldPrefix;
 631             }
 632         }
 633     }
 634 }
 635
 636
 637 void Foam::refinementHistory::distribute(const mapDistributePolyMesh& map)
 638 {
 639     if (!active())
 640     {
 641         FatalErrorIn
 642         (
 643             "refinementHistory::distribute(const mapDistributePolyMesh&)"
 644         )   << "Calling distribute on inactive history" << abort(FatalError);
 645     }
 646
 647
 648     if (!Pstream::parRun())
 649     {
 650         return;
 651     }
 652
 653     // Remove unreferenced history.
 654     compact();
 655
 656     Pout<< nl << "--BEFORE:" << endl;
 657     writeDebug();
 658     Pout<< "---------" << nl << endl;
 659
 660
 661     // Distribution is only partially functional.
 662     // If all 8 cells resulting from a single parent are sent across in one
 663     // go it will also send across that part of the refinement history.
 664     // If however e.g. first 1 and then the other 7 are sent across the
 665     // history will not be reconstructed.
 666
 667     // Determine clusters. This is per every entry in splitCells_ (that is
 668     // a parent of some refinement) a label giving the processor it goes to
 669     // if all its children are going to the same processor.
 670
 671     // Per visible cell the processor it goes to.
 672     labelList destination(visibleCells_.size());
 673
 674     const labelListList& subCellMap = map.cellMap().subMap();
 675
 676     forAll(subCellMap, procI)
 677     {
 678         const labelList& newToOld = subCellMap[procI];
 679
 680         forAll(newToOld, i)
 681         {
 682             label oldCellI = newToOld[i];
 683
 684             destination[oldCellI] = procI;
 685         }
 686     }
 687
 688 //Pout<< "refinementHistory::distribute :"
 689 //    << " destination:" << destination << endl;
 690
 691     // Per splitCell entry the processor it moves to
 692     labelList splitCellProc(splitCells_.size(), -1);
 693     // Per splitCell entry the number of live cells that move to that processor
 694     labelList splitCellNum(splitCells_.size(), 0);
 695
 696     forAll(visibleCells_, cellI)
 697     {
 698         label index = visibleCells_[cellI];
 699
 700         if (index >= 0)
 701         {
 702             countProc
 703             (
 704                 splitCells_[index].parent_,
 705                 destination[cellI],
 706                 splitCellProc,
 707                 splitCellNum
 708             );
 709         }
 710     }
 711
 712 Pout<< "refinementHistory::distribute :"
 713     << " splitCellProc:" << splitCellProc << endl;
 714
 715 Pout<< "refinementHistory::distribute :"
 716     << " splitCellNum:" << splitCellNum << endl;
 717
 718
 719     // Create subsetted refinement tree consisting of all parents that
 720     // move in their whole to other processor.
 721     for (label procI = 0; procI < Pstream::nProcs(); procI++)
 722     {
 723         Pout<< "-- Subetting for processor " << procI << endl;
 724
 725         // From uncompacted to compacted splitCells.
 726         labelList oldToNew(splitCells_.size(), -1);
 727
 728         // Compacted splitCells. Similar to subset routine below.
 729         DynamicList<splitCell8> newSplitCells(splitCells_.size());
 730
 731         // Loop over all entries. Note: could recurse like countProc so only
 732         // visit used entries but is probably not worth it.
 733
 734         forAll(splitCells_, index)
 735         {
 736 //            Pout<< "oldCell:" << index
 737 //                << " proc:" << splitCellProc[index]
 738 //                << " nCells:" << splitCellNum[index]
 739 //                << endl;
 740
 741             if (splitCellProc[index] == procI && splitCellNum[index] == 8)
 742             {
 743                 // Entry moves in its whole to procI
 744                 oldToNew[index] = newSplitCells.size();
 745                 newSplitCells.append(splitCells_[index]);
 746
 747                 Pout<< "Added oldCell " << index
 748                     << " info " << newSplitCells[newSplitCells.size()-1]
 749                     << " at position " << newSplitCells.size()-1
 750                     << endl;
 751             }
 752         }
 753
 754         // Add live cells that are subsetted.
 755         forAll(visibleCells_, cellI)
 756         {
 757             label index = visibleCells_[cellI];
 758
 759             if (index >= 0 && destination[cellI] == procI)
 760             {
 761                 label parent = splitCells_[index].parent_;
 762
 763                 Pout<< "Adding refined cell " << cellI
 764                     << " since moves to "
 765                     << procI << " old parent:" << parent << endl;
 766
 767                 // Create new splitCell with parent
 768                 oldToNew[index] = newSplitCells.size();
 769                 newSplitCells.append(splitCell8(parent));
 770             }
 771         }
 772
 773         //forAll(oldToNew, index)
 774         //{
 775         //    Pout<< "old:" << index << " new:" << oldToNew[index]
 776         //        << endl;
 777         //}
 778
 779         newSplitCells.shrink();
 780
 781         // Renumber contents of newSplitCells
 782         forAll(newSplitCells, index)
 783         {
 784             splitCell8& split = newSplitCells[index];
 785
 786             if (split.parent_ >= 0)
 787             {
 788                 split.parent_ = oldToNew[split.parent_];
 789             }
 790             if (split.addedCellsPtr_.valid())
 791             {
 792                 FixedList<label, 8>& splits = split.addedCellsPtr_();
 793
 794                 forAll(splits, i)
 795                 {
 796                     if (splits[i] >= 0)
 797                     {
 798                         splits[i] = oldToNew[splits[i]];
 799                     }
 800                 }
 801             }
 802         }
 803
 804
 805         const labelList& subMap = subCellMap[procI];
 806
 807         // New visible cells.
 808         labelList newVisibleCells(subMap.size(), -1);
 809
 810         forAll(subMap, newCellI)
 811         {
 812             label oldCellI = subMap[newCellI];
 813
 814             label oldIndex = visibleCells_[oldCellI];
 815
 816             if (oldIndex >= 0)
 817             {
 818                 newVisibleCells[newCellI] = oldToNew[oldIndex];
 819             }
 820         }
 821
 822         //Pout<< nl << "--Subset for domain:" << procI << endl;
 823         //writeDebug(newVisibleCells, newSplitCells);
 824         //Pout<< "---------" << nl << endl;
 825
 826
 827         // Send to neighbours
 828         OPstream toNbr(Pstream::blocking, procI);
 829         toNbr << newSplitCells << newVisibleCells;
 830     }
 831
 832
 833     // Receive from neighbours and merge
 834     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 835
 836     // Remove all entries. Leave storage intact.
 837     splitCells_.clear();
 838
 839     visibleCells_.setSize(map.mesh().nCells());
 840     visibleCells_ = -1;
 841
 842     for (label procI = 0; procI < Pstream::nProcs(); procI++)
 843     {
 844         IPstream fromNbr(Pstream::blocking, procI);
 845         List<splitCell8> newSplitCells(fromNbr);
 846         labelList newVisibleCells(fromNbr);
 847
 848         //Pout<< nl << "--Received from domain:" << procI << endl;
 849         //writeDebug(newVisibleCells, newSplitCells);
 850         //Pout<< "---------" << nl << endl;
 851
 852
 853         // newSplitCells contain indices only into newSplitCells so
 854         // renumbering can be done here.
 855         label offset = splitCells_.size();
 856
 857         Pout<< "**Renumbering data from proc " << procI << " with offset "
 858             << offset << endl;
 859
 860         forAll(newSplitCells, index)
 861         {
 862             splitCell8& split = newSplitCells[index];
 863
 864             if (split.parent_ >= 0)
 865             {
 866                 split.parent_ += offset;
 867             }
 868             if (split.addedCellsPtr_.valid())
 869             {
 870                 FixedList<label, 8>& splits = split.addedCellsPtr_();
 871
 872                 forAll(splits, i)
 873                 {
 874                     if (splits[i] >= 0)
 875                     {
 876                         splits[i] += offset;
 877                     }
 878                 }
 879             }
 880
 881             splitCells_.append(split);
 882         }
 883
 884
 885         // Combine visibleCell.
 886         const labelList& constructMap = map.cellMap().constructMap()[procI];
 887
 888         forAll(newVisibleCells, i)
 889         {
 890             visibleCells_[constructMap[i]] = newVisibleCells[i] + offset;
 891         }
 892     }
 893     splitCells_.shrink();
 894
 895     Pout<< nl << "--AFTER:" << endl;
 896     writeDebug();
 897     Pout<< "---------" << nl << endl;
 898 }
 899
 900
 901 void Foam::refinementHistory::compact()
 902 {
 903     if (debug)
 904     {
 905         Pout<< "refinementHistory::compact() Entering with:"
 906             << " freeSplitCells_:" << freeSplitCells_.size()
 907             << " splitCells_:" << splitCells_.size()
 908             << " visibleCells_:" << visibleCells_.size()
 909             << endl;
 910
 911         // Check all free splitCells are marked as such
 912         forAll(freeSplitCells_, i)
 913         {
 914             label index = freeSplitCells_[i];
 915
 916             if (splitCells_[index].parent_ != -2)
 917             {
 918                 FatalErrorIn("refinementHistory::compact()")
 919                     << "Problem index:" << index
 920                     << abort(FatalError);
 921             }
 922         }
 923
 924         // Check none of the visible cells are marked as free
 925         forAll(visibleCells_, cellI)
 926         {
 927             if
 928             (
 929                 visibleCells_[cellI] >= 0
 930              && splitCells_[visibleCells_[cellI]].parent_ == -2
 931             )
 932             {
 933                 FatalErrorIn("refinementHistory::compact()")
 934                     << "Problem : visible cell:" << cellI
 935                     << " is marked as being free." << abort(FatalError);
 936             }
 937         }
 938     }
 939
 940     DynamicList<splitCell8> newSplitCells(splitCells_.size());
 941
 942     // From uncompacted to compacted splitCells.
 943     labelList oldToNew(splitCells_.size(), -1);
 944
 945     // Mark all used splitCell entries. These are either indexed by visibleCells
 946     // or indexed from other splitCell entries.
 947
 948     // Mark from visibleCells
 949     forAll(visibleCells_, cellI)
 950     {
 951         label index = visibleCells_[cellI];
 952
 953         if (index >= 0)
 954         {
 955             // Make sure we only mark visible indices if they either have a
 956             // parent or subsplits.
 957             if
 958             (
 959                 splitCells_[index].parent_ != -1
 960              || splitCells_[index].addedCellsPtr_.valid()
 961             )
 962             {
 963                 markSplit(index, oldToNew, newSplitCells);
 964             }
 965         }
 966     }
 967
 968     // Mark from splitCells
 969     forAll(splitCells_, index)
 970     {
 971         if (splitCells_[index].parent_ == -2)
 972         {
 973             // freed cell.
 974         }
 975         else if
 976         (
 977             splitCells_[index].parent_ == -1
 978          && splitCells_[index].addedCellsPtr_.empty()
 979         )
 980         {
 981             // recombined cell. No need to keep since no parent and no subsplits
 982             // Note that gets marked if reachable from other index!
 983         }
 984         else
 985         {
 986             // Is used element.
 987             markSplit(index, oldToNew, newSplitCells);
 988         }
 989     }
 990
 991
 992     // Now oldToNew is fully complete and compacted elements are in
 993     // newSplitCells.
 994     // Renumber contents of newSplitCells and visibleCells.
 995     forAll(newSplitCells, index)
 996     {
 997         splitCell8& split = newSplitCells[index];
 998
 999         if (split.parent_ >= 0)
1000         {
1001             split.parent_ = oldToNew[split.parent_];
1002         }
1003         if (split.addedCellsPtr_.valid())
1004         {
1005             FixedList<label, 8>& splits = split.addedCellsPtr_();
1006
1007             forAll(splits, i)
1008             {
1009                 if (splits[i] >= 0)
1010                 {
1011                     splits[i] = oldToNew[splits[i]];
1012                 }
1013             }
1014         }
1015     }
1016
1017
1018     if (debug)
1019     {
1020         Pout<< "refinementHistory::compact : compacted splitCells from "
1021             << splitCells_.size() << " to " << newSplitCells.size() << endl;
1022     }
1023
1024     splitCells_.transfer(newSplitCells);
1025     freeSplitCells_.clearStorage();
1026
1027
1028     if (debug)
1029     {
1030         Pout<< "refinementHistory::compact() NOW:"
1031             << " freeSplitCells_:" << freeSplitCells_.size()
1032             << " splitCells_:" << splitCells_.size()
1033             << " newSplitCells:" << newSplitCells.size()
1034             << " visibleCells_:" << visibleCells_.size()
1035             << endl;
1036     }
1037
1038
1039     // Adapt indices in visibleCells_
1040     forAll(visibleCells_, cellI)
1041     {
1042         label index = visibleCells_[cellI];
1043
1044         if (index >= 0)
1045         {
1046             // Note that oldToNew can be -1 so it resets newVisibleCells.
1047             visibleCells_[cellI] = oldToNew[index];
1048         }
1049         else
1050         {
1051             // Keep -1 value.
1052         }
1053     }
1054 }
1055
1056
1057 void Foam::refinementHistory::writeDebug() const
1058 {
1059     writeDebug(visibleCells_, splitCells_);
1060 }
1061
1062
1063 void Foam::refinementHistory::storeSplit
1064 (
1065     const label cellI,
1066     const labelList& addedCells
1067 )
1068 {
1069     label parentIndex = -1;
1070
1071     if (visibleCells_[cellI] != -1)
1072     {
1073         // Was already live. The current live cell becomes the
1074         // parent of the cells split off from it.
1075
1076         parentIndex = visibleCells_[cellI];
1077
1078         // It is no longer live (note that actually cellI gets alive
1079         // again below since is addedCells[0])
1080         visibleCells_[cellI] = -1;
1081     }
1082     else
1083     {
1084         // Create 0th level. -1 parent to denote this.
1085         parentIndex = allocateSplitCell(-1, -1);
1086     }
1087
1088     // Create live entries for added cells that point to the
1089     // cell they were created from (parentIndex)
1090     forAll(addedCells, i)
1091     {
1092         label addedCellI = addedCells[i];
1093
1094         // Create entries for the split off cells. All of them
1095         // are visible.
1096         visibleCells_[addedCellI] = allocateSplitCell(parentIndex, i);
1097     }
1098 }
1099
1100
1101 void Foam::refinementHistory::combineCells
1102 (
1103     const label masterCellI,
1104     const labelList& combinedCells
1105 )
1106 {
1107     // Save the parent structure
1108     label parentIndex = splitCells_[visibleCells_[masterCellI]].parent_;
1109
1110     // Remove the information for the combined cells
1111     forAll(combinedCells, i)
1112     {
1113         label cellI = combinedCells[i];
1114
1115         freeSplitCell(visibleCells_[cellI]);
1116         visibleCells_[cellI] = -1;
1117     }
1118
1119     splitCell8& parentSplit = splitCells_[parentIndex];
1120     parentSplit.addedCellsPtr_.reset(NULL);
1121     visibleCells_[masterCellI] = parentIndex;
1122 }
1123
1124
1125 bool Foam::refinementHistory::readData(Istream& is)
1126 {
1127     is >> *this;
1128     return !is.bad();
1129 }
1130
1131
1132 bool Foam::refinementHistory::writeData(Ostream& os) const
1133 {
1134     os << *this;
1135
1136     return os.good();
1137 }
1138
1139
1140 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
1141
1142 Foam::Istream& Foam::operator>>(Istream& is, refinementHistory& rh)
1143 {
1144     rh.freeSplitCells_.clearStorage();
1145
1146     is >> rh.splitCells_ >> rh.visibleCells_;
1147
1148     // Check indices.
1149     rh.checkIndices();
1150
1151     return is;
1152 }
1153
1154
1155 Foam::Ostream& Foam::operator<<(Ostream& os, const refinementHistory& rh)
1156 {
1157     const_cast<refinementHistory&>(rh).compact();
1158
1159     return os   << "// splitCells" << nl
1160                 << rh.splitCells_ << nl
1161                 << "// visibleCells" << nl
1162                 << rh.visibleCells_;
1163 }
1164
1165
1166 // ************************************************************************* //