src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/dynamicTopoFvMesh.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   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 Class
  26     dynamicTopoFvMesh
  27
  28 Description
  29     An implementation of dynamic changes to mesh-topology
  30
  31 Author
  32     Sandeep Menon
  33     University of Massachusetts Amherst
  34     All rights reserved
  35
  36 \*---------------------------------------------------------------------------*/
  37
  38 #include "dynamicTopoFvMesh.H"
  39 #include "addToRunTimeSelectionTable.H"
  40
  41 #include "eMesh.H"
  42 #include "Stack.H"
  43 #include "triFace.H"
  44 #include "changeMap.H"
  45 #include "clockTime.H"
  46 #include "volFields.H"
  47 #include "MeshObject.H"
  48 #include "topoMapper.H"
  49 #include "coupledInfo.H"
  50 #include "mapPolyMesh.H"
  51 #include "MapFvFields.H"
  52 #include "SortableList.H"
  53 #include "motionSolver.H"
  54 #include "fvPatchFields.H"
  55 #include "fvsPatchFields.H"
  56 #include "lengthScaleEstimator.H"
  57 #include "conservativeMapFields.H"
  58
  59 namespace Foam
  60 {
  61
  62 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  63
  64 defineTypeNameAndDebug(dynamicTopoFvMesh,0);
  65
  66 addToRunTimeSelectionTable(dynamicFvMesh, dynamicTopoFvMesh, IOobject);
  67
  68 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  69
  70 // Construct from IOobject
  71 dynamicTopoFvMesh::dynamicTopoFvMesh(const IOobject& io)
  72 :
  73     dynamicFvMesh(io),
  74     baseMesh_(*this),
  75     nPatches_(boundaryMesh().size()),
  76     topoChangeFlag_(false),
  77     isSubMesh_(false),
  78     dict_
  79     (
  80         IOobject
  81         (
  82             "dynamicMeshDict",
  83             polyMesh::time().constant(),
  84             (*this),
  85             IOobject::MUST_READ,
  86             IOobject::NO_WRITE,
  87             false
  88         )
  89     ),
  90     mandatory_
  91     (
  92         dict_.subDict("dynamicTopoFvMesh").lookup("allOptionsMandatory")
  93     ),
  94     twoDMesh_(polyMesh::nGeometricD() == 2 ? true : false),
  95     edgeRefinement_
  96     (
  97         dict_.subDict("dynamicTopoFvMesh").lookup("edgeRefinement")
  98     ),
  99     loadMotionSolver_(true),
 100     bandWidthReduction_(false),
 101     coupledModification_(false),
 102     interval_(1),
 103     eMeshPtr_(NULL),
 104     mapper_(NULL),
 105     motionSolver_(NULL),
 106     lengthEstimator_(NULL),
 107     oldPoints_(polyMesh::points()),
 108     points_(polyMesh::points()),
 109     faces_(polyMesh::faces()),
 110     owner_(polyMesh::faceOwner()),
 111     neighbour_(polyMesh::faceNeighbour()),
 112     cells_(primitiveMesh::cells()),
 113     oldPatchSizes_(nPatches_, 0),
 114     patchSizes_(nPatches_, 0),
 115     oldPatchStarts_(nPatches_, -1),
 116     patchStarts_(nPatches_, -1),
 117     oldEdgePatchSizes_(nPatches_, 0),
 118     edgePatchSizes_(nPatches_, 0),
 119     oldEdgePatchStarts_(nPatches_, -1),
 120     edgePatchStarts_(nPatches_, -1),
 121     oldPatchNMeshPoints_(nPatches_, -1),
 122     patchNMeshPoints_(nPatches_, -1),
 123     nOldPoints_(primitiveMesh::nPoints()),
 124     nPoints_(primitiveMesh::nPoints()),
 125     nOldEdges_(0),
 126     nEdges_(0),
 127     nOldFaces_(primitiveMesh::nFaces()),
 128     nFaces_(primitiveMesh::nFaces()),
 129     nOldCells_(primitiveMesh::nCells()),
 130     nCells_(primitiveMesh::nCells()),
 131     nOldInternalFaces_(primitiveMesh::nInternalFaces()),
 132     nInternalFaces_(primitiveMesh::nInternalFaces()),
 133     nOldInternalEdges_(0),
 134     nInternalEdges_(0),
 135     maxModifications_(-1),
 136     statistics_(0),
 137     sliverThreshold_(0.1),
 138     slicePairs_(0),
 139     maxTetsPerEdge_(-1),
 140     swapDeviation_(0.0),
 141     allowTableResize_(false)
 142 {
 143     // Check the size of owner/neighbour
 144     if (owner_.size() != neighbour_.size())
 145     {
 146         // Size up to number of faces
 147         neighbour_.setSize(nFaces_, -1);
 148     }
 149
 150     const polyBoundaryMesh& boundary = polyMesh::boundaryMesh();
 151
 152     // Initialize the multiThreading environment
 153     initializeThreadingEnvironment();
 154
 155     // Read optional parameters.
 156     readOptionalParameters();
 157
 158     // Initialize patch-size information
 159     for (label i = 0; i < nPatches_; i++)
 160     {
 161         patchNMeshPoints_[i] = boundary[i].meshPoints().size();
 162         oldPatchSizes_[i] = patchSizes_[i] = boundary[i].size();
 163         oldPatchStarts_[i] = patchStarts_[i] = boundary[i].start();
 164         oldPatchNMeshPoints_[i] = patchNMeshPoints_[i];
 165     }
 166
 167     // Open the tetMetric dynamic-link library (for 3D only)
 168     loadMetric();
 169
 170     // Initialize edge-related connectivity structures
 171     initEdges();
 172
 173     // Load the mesh-motion solver
 174     loadMotionSolver();
 175
 176     // Load the field-mapper
 177     loadFieldMapper();
 178
 179     // Set sizes for the reverse maps
 180     reversePointMap_.setSize(nPoints_, -7);
 181     reverseEdgeMap_.setSize(nEdges_, -7);
 182     reverseFaceMap_.setSize(nFaces_, -7);
 183     reverseCellMap_.setSize(nCells_, -7);
 184
 185     // Load the length-scale estimator,
 186     // and read refinement options
 187     loadLengthScaleEstimator();
 188 }
 189
 190
 191 //- Construct from components. Used for subMeshes only.
 192 dynamicTopoFvMesh::dynamicTopoFvMesh
 193 (
 194     const dynamicTopoFvMesh& mesh,
 195     const IOobject& io,
 196     const Xfer<pointField>& points,
 197     const Xfer<pointField>& oldPoints,
 198     const Xfer<edgeList>& edges,
 199     const Xfer<faceList>& faces,
 200     const Xfer<labelListList>& faceEdges,
 201     const Xfer<labelList>& owner,
 202     const Xfer<labelList>& neighbour,
 203     const labelList& faceStarts,
 204     const labelList& faceSizes,
 205     const labelList& edgeStarts,
 206     const labelList& edgeSizes,
 207     const wordList& patchNames,
 208     const dictionary& patchDict
 209 )
 210 :
 211     dynamicFvMesh
 212     (
 213         io,
 214         oldPoints,
 215         faces,
 216         owner,
 217         neighbour,
 218         false
 219     ),
 220     baseMesh_(mesh),
 221     nPatches_(faceStarts.size()),
 222     topoChangeFlag_(false),
 223     isSubMesh_(true),
 224     dict_(mesh.dict_),
 225     mandatory_(mesh.mandatory_),
 226     twoDMesh_(mesh.twoDMesh_),
 227     edgeRefinement_(mesh.edgeRefinement_),
 228     loadMotionSolver_(mesh.loadMotionSolver_),
 229     bandWidthReduction_(mesh.bandWidthReduction_),
 230     coupledModification_(false),
 231     interval_(1),
 232     eMeshPtr_(NULL),
 233     mapper_(NULL),
 234     motionSolver_(NULL),
 235     lengthEstimator_(NULL),
 236     oldPoints_(polyMesh::points()),
 237     points_(points()),
 238     faces_(polyMesh::faces()),
 239     cells_(polyMesh::cells()),
 240     edges_(edges),
 241     faceEdges_(faceEdges),
 242     oldPatchSizes_(faceSizes),
 243     patchSizes_(faceSizes),
 244     oldPatchStarts_(faceStarts),
 245     patchStarts_(faceStarts),
 246     oldEdgePatchSizes_(edgeSizes),
 247     edgePatchSizes_(edgeSizes),
 248     oldEdgePatchStarts_(edgeStarts),
 249     edgePatchStarts_(edgeStarts),
 250     oldPatchNMeshPoints_(faceStarts.size(), -1),
 251     patchNMeshPoints_(faceStarts.size(), -1),
 252     nOldPoints_(points_.size()),
 253     nPoints_(points_.size()),
 254     nOldEdges_(edges_.size()),
 255     nEdges_(edges_.size()),
 256     nOldFaces_(faces_.size()),
 257     nFaces_(faces_.size()),
 258     nOldCells_(cells_.size()),
 259     nCells_(cells_.size()),
 260     nOldInternalFaces_(faceStarts[0]),
 261     nInternalFaces_(faceStarts[0]),
 262     nOldInternalEdges_(edgeStarts[0]),
 263     nInternalEdges_(edgeStarts[0]),
 264     maxModifications_(mesh.maxModifications_),
 265     statistics_(0),
 266     sliverThreshold_(mesh.sliverThreshold_),
 267     slicePairs_(0),
 268     maxTetsPerEdge_(mesh.maxTetsPerEdge_),
 269     swapDeviation_(mesh.swapDeviation_),
 270     allowTableResize_(mesh.allowTableResize_),
 271     tetMetric_(mesh.tetMetric_)
 272 {
 273     // Initialize owner and neighbour
 274     owner_.setSize(faces_.size(), -1);
 275     neighbour_.setSize(faces_.size(), -1);
 276
 277     // Set owner and neighbour from polyMesh
 278     const labelList& own = polyMesh::faceOwner();
 279     const labelList& nei = polyMesh::faceNeighbour();
 280
 281     forAll(own, faceI)
 282     {
 283         owner_[faceI] = own[faceI];
 284     }
 285
 286     forAll(nei, faceI)
 287     {
 288         neighbour_[faceI] = nei[faceI];
 289     }
 290
 291     // Initialize the multiThreading environment.
 292     // Force to single-threaded.
 293     initializeThreadingEnvironment(1);
 294
 295     const polyBoundaryMesh& boundary = polyMesh::boundaryMesh();
 296
 297     // Add a boundary patches for polyMesh.
 298     List<polyPatch*> patches(faceStarts.size());
 299
 300     forAll(patches, patchI)
 301     {
 302         // Set the pointer
 303         patches[patchI] =
 304         (
 305             polyPatch::New
 306             (
 307                 patchNames[patchI],
 308                 patchDict.subDict
 309                 (
 310                     patchNames[patchI]
 311                 ),
 312                 patchI,
 313                 boundary
 314             ).ptr()
 315         );
 316     }
 317
 318     // Add patches, but don't calculate geometry, etc
 319     fvMesh::addFvPatches(patches, false);
 320
 321     // Set sizes for the reverse maps
 322     reversePointMap_.setSize(nPoints_, -7);
 323     reverseEdgeMap_.setSize(nEdges_, -7);
 324     reverseFaceMap_.setSize(nFaces_, -7);
 325     reverseCellMap_.setSize(nCells_, -7);
 326
 327     // Now build edgeFaces and pointEdges information.
 328     edgeFaces_ = invertManyToMany<labelList, labelList>(nEdges_, faceEdges_);
 329
 330     if (!twoDMesh_)
 331     {
 332         pointEdges_ = invertManyToMany<edge, labelList>(nPoints_, edges_);
 333     }
 334 }
 335
 336
 337 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 338
 339 dynamicTopoFvMesh::~dynamicTopoFvMesh()
 340 {}
 341
 342
 343 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 344
 345 // Insert the specified cell to the mesh.
 346 label dynamicTopoFvMesh::insertCell
 347 (
 348     const cell& newCell,
 349     const scalar lengthScale,
 350     const label zoneID
 351 )
 352 {
 353     label newCellIndex = cells_.size();
 354
 355     if (debug > 2)
 356     {
 357         Pout<< "Inserting cell: "
 358             << newCellIndex << ": "
 359             << newCell << endl;
 360     }
 361
 362     cells_.append(newCell);
 363
 364     if (edgeRefinement_)
 365     {
 366         lengthScale_.append(lengthScale);
 367     }
 368
 369     // Add to the zone if necessary
 370     if (zoneID >= 0)
 371     {
 372         addedCellZones_.insert(newCellIndex, zoneID);
 373     }
 374
 375     nCells_++;
 376
 377     return newCellIndex;
 378 }
 379
 380
 381 // Remove the specified cell from the mesh
 382 void dynamicTopoFvMesh::removeCell
 383 (
 384     const label cIndex
 385 )
 386 {
 387     if (debug > 2)
 388     {
 389         Pout<< "Removing cell: "
 390             << cIndex << ": "
 391             << cells_[cIndex]
 392             << endl;
 393     }
 394
 395     cells_[cIndex].clear();
 396
 397     if (edgeRefinement_)
 398     {
 399         lengthScale_[cIndex] = -1.0;
 400     }
 401
 402     // Update the number of cells, and the reverse cell map
 403     nCells_--;
 404
 405     if (cIndex < nOldCells_)
 406     {
 407         reverseCellMap_[cIndex] = -1;
 408     }
 409     else
 410     {
 411         // Store this information for the reOrdering stage
 412         deletedCells_.insert(cIndex);
 413     }
 414
 415     // Check if this cell was added to a zone
 416     Map<label>::iterator it = addedCellZones_.find(cIndex);
 417
 418     if (it != addedCellZones_.end())
 419     {
 420         addedCellZones_.erase(it);
 421     }
 422
 423     // Check if the cell was added in the current morph, and delete
 424     forAll(cellsFromPoints_, indexI)
 425     {
 426         if (cellsFromPoints_[indexI].index() == cIndex)
 427         {
 428             // Remove entry from the list
 429             meshOps::removeIndex(indexI, cellsFromPoints_);
 430             break;
 431         }
 432     }
 433
 434     forAll(cellsFromEdges_, indexI)
 435     {
 436         if (cellsFromEdges_[indexI].index() == cIndex)
 437         {
 438             // Remove entry from the list
 439             meshOps::removeIndex(indexI, cellsFromEdges_);
 440             break;
 441         }
 442     }
 443
 444     forAll(cellsFromFaces_, indexI)
 445     {
 446         if (cellsFromFaces_[indexI].index() == cIndex)
 447         {
 448             // Remove entry from the list
 449             meshOps::removeIndex(indexI, cellsFromFaces_);
 450             break;
 451         }
 452     }
 453
 454     forAll(cellsFromCells_, indexI)
 455     {
 456         if (cellsFromCells_[indexI].index() == cIndex)
 457         {
 458             // Remove entry from the list
 459             meshOps::removeIndex(indexI, cellsFromCells_);
 460             break;
 461         }
 462     }
 463 }
 464
 465
 466 // Utility method for face-insertion
 467 label dynamicTopoFvMesh::insertFace
 468 (
 469     const label patch,
 470     const face& newFace,
 471     const label newOwner,
 472     const label newNeighbour,
 473     const label zoneID
 474 )
 475 {
 476     // Append the specified face to each face-related list.
 477     // Reordering is performed after all pending changes
 478     // (flips, bisections, contractions, etc) have been made to the mesh
 479     label newFaceIndex = faces_.size();
 480
 481     faces_.append(newFace);
 482     owner_.append(newOwner);
 483     neighbour_.append(newNeighbour);
 484
 485     if (debug > 2)
 486     {
 487         Pout<< "Inserting face: "
 488             << newFaceIndex << ": "
 489             << newFace
 490             << " Owner: " << newOwner
 491             << " Neighbour: " << newNeighbour;
 492
 493         const polyBoundaryMesh& boundary = boundaryMesh();
 494
 495         Pout<< " Patch: ";
 496
 497         if (patch == -1)
 498         {
 499             Pout<< "Internal" << endl;
 500         }
 501         else
 502         if (patch < boundary.size())
 503         {
 504             Pout<< boundary[patch].name() << endl;
 505         }
 506         else
 507         {
 508             Pout<< " New patch: " << patch << endl;
 509         }
 510     }
 511
 512     // Keep track of added boundary faces in a separate hash-table
 513     // This information will be required at the reordering stage
 514     addedFacePatches_.insert(newFaceIndex,patch);
 515
 516     if (newNeighbour == -1)
 517     {
 518         // Modify patch information for this boundary face
 519         patchSizes_[patch]++;
 520
 521         for (label i = (patch + 1); i < nPatches_; i++)
 522         {
 523             patchStarts_[i]++;
 524         }
 525     }
 526     else
 527     {
 528         // Increment the number of internal faces,
 529         // and subsequent patch-starts
 530         nInternalFaces_++;
 531
 532         for (label i = 0; i < nPatches_; i++)
 533         {
 534             patchStarts_[i]++;
 535         }
 536     }
 537
 538     // Add to the zone if explicitly specified
 539     if (zoneID >= 0)
 540     {
 541         addedFaceZones_.insert(newFaceIndex, zoneID);
 542     }
 543     else
 544     {
 545         // No zone was specified. Check if this
 546         // face is added to a coupled patch associated
 547         // with a faceZone.
 548         forAll(patchCoupling_, patchI)
 549         {
 550             if (patchCoupling_(patchI))
 551             {
 552                 if (patchI == patch)
 553                 {
 554                     if (patchCoupling_[patchI].masterFaceZone() > -1)
 555                     {
 556                         addedFaceZones_.insert
 557                         (
 558                             newFaceIndex,
 559                             patchCoupling_[patchI].masterFaceZone()
 560                         );
 561                     }
 562
 563                     break;
 564                 }
 565
 566                 if (patchCoupling_[patchI].map().slaveIndex() == patch)
 567                 {
 568                     if (patchCoupling_[patchI].slaveFaceZone() > -1)
 569                     {
 570                         addedFaceZones_.insert
 571                         (
 572                             newFaceIndex,
 573                             patchCoupling_[patchI].slaveFaceZone()
 574                         );
 575                     }
 576
 577                     break;
 578                 }
 579             }
 580         }
 581     }
 582
 583     // Increment the total face count
 584     nFaces_++;
 585
 586     return newFaceIndex;
 587 }
 588
 589
 590 // Remove the specified face from the mesh
 591 void dynamicTopoFvMesh::removeFace
 592 (
 593     const label fIndex
 594 )
 595 {
 596     // Identify the patch for this face
 597     label patch = whichPatch(fIndex);
 598
 599     if (debug > 2)
 600     {
 601         Pout<< "Removing face: "
 602             << fIndex << ": "
 603             << faces_[fIndex];
 604
 605         const polyBoundaryMesh& boundary = boundaryMesh();
 606
 607         Pout<< " Patch: ";
 608
 609         if (patch == -1)
 610         {
 611             Pout<< "Internal" << endl;
 612         }
 613         else
 614         if (patch < boundary.size())
 615         {
 616             Pout<< boundary[patch].name() << endl;
 617         }
 618         else
 619         {
 620             Pout<< " New patch: " << patch << endl;
 621         }
 622     }
 623
 624     if (patch >= 0)
 625     {
 626         // Modify patch information for this boundary face
 627         patchSizes_[patch]--;
 628
 629         for (label i = (patch + 1); i < nPatches_; i++)
 630         {
 631             patchStarts_[i]--;
 632         }
 633     }
 634     else
 635     {
 636         // Decrement the internal face count, and subsequent patch-starts
 637         nInternalFaces_--;
 638
 639         forAll(patchStarts_, patchI)
 640         {
 641             patchStarts_[patchI]--;
 642         }
 643     }
 644
 645     // Clear entities.
 646     faces_[fIndex].clear();
 647     owner_[fIndex] = -1;
 648     neighbour_[fIndex] = -1;
 649     faceEdges_[fIndex].clear();
 650
 651     // Entity won't be removed from the stack for efficiency
 652     // It will be discarded on access instead.
 653
 654     // Update coupled face maps, if necessary.
 655     forAll(patchCoupling_, patchI)
 656     {
 657         if (patchCoupling_(patchI))
 658         {
 659             const coupleMap& cMap = patchCoupling_[patchI].map();
 660
 661             cMap.removeMaster(coupleMap::FACE, fIndex);
 662             cMap.removeSlave(coupleMap::FACE, fIndex);
 663         }
 664     }
 665
 666     // Update the reverse face-map, but only if this is a face that existed
 667     // at time [n]. Added faces which are deleted during the topology change
 668     // needn't be updated.
 669     if (fIndex < nOldFaces_)
 670     {
 671         reverseFaceMap_[fIndex] = -1;
 672     }
 673     else
 674     {
 675         // Store this information for the reOrdering stage
 676         deletedFaces_.insert(fIndex);
 677     }
 678
 679     // Check and remove from the list of added face patches
 680     Map<label>::iterator fpit = addedFacePatches_.find(fIndex);
 681
 682     if (fpit != addedFacePatches_.end())
 683     {
 684         addedFacePatches_.erase(fpit);
 685     }
 686
 687     // Check if this face was added to a zone
 688     Map<label>::iterator fzit = addedFaceZones_.find(fIndex);
 689
 690     if (fzit != addedFaceZones_.end())
 691     {
 692         addedFaceZones_.erase(fzit);
 693     }
 694
 695     // Check if the face was added in the current morph, and delete
 696     forAll(facesFromPoints_, indexI)
 697     {
 698         if (facesFromPoints_[indexI].index() == fIndex)
 699         {
 700             // Remove entry from the list
 701             meshOps::removeIndex
 702             (
 703                 indexI,
 704                 facesFromPoints_
 705             );
 706
 707             break;
 708         }
 709     }
 710
 711     forAll(facesFromEdges_, indexI)
 712     {
 713         if (facesFromEdges_[indexI].index() == fIndex)
 714         {
 715             // Remove entry from the list
 716             meshOps::removeIndex
 717             (
 718                 indexI,
 719                 facesFromEdges_
 720             );
 721
 722             break;
 723         }
 724     }
 725
 726     forAll(facesFromFaces_, indexI)
 727     {
 728         if (facesFromFaces_[indexI].index() == fIndex)
 729         {
 730             // Remove entry from the list
 731             meshOps::removeIndex
 732             (
 733                 indexI,
 734                 facesFromFaces_
 735             );
 736
 737             break;
 738         }
 739     }
 740
 741     // Remove from the flipFaces list, if necessary
 742     labelHashSet::iterator ffit = flipFaces_.find(fIndex);
 743
 744     if (ffit != flipFaces_.end())
 745     {
 746         flipFaces_.erase(ffit);
 747     }
 748
 749     // Decrement the total face-count
 750     nFaces_--;
 751 }
 752
 753
 754 // Insert the specified edge to the mesh
 755 label dynamicTopoFvMesh::insertEdge
 756 (
 757     const label patch,
 758     const edge& newEdge,
 759     const labelList& edgeFaces
 760 )
 761 {
 762     label newEdgeIndex = edges_.size();
 763
 764     edges_.append(newEdge);
 765     edgeFaces_.append(edgeFaces);
 766
 767     if (debug > 2)
 768     {
 769         Pout<< "Inserting edge: "
 770             << newEdgeIndex << ": "
 771             << newEdge;
 772
 773         const polyBoundaryMesh& boundary = boundaryMesh();
 774
 775         Pout<< " Patch: ";
 776
 777         if (patch == -1)
 778         {
 779             Pout<< "Internal" << endl;
 780         }
 781         else
 782         if (patch < boundary.size())
 783         {
 784             Pout<< boundary[patch].name() << endl;
 785         }
 786         else
 787         {
 788             Pout<< " New patch: " << patch << endl;
 789         }
 790     }
 791
 792     // Keep track of added edges in a separate hash-table
 793     // This information will be required at the reordering stage
 794     addedEdgePatches_.insert(newEdgeIndex,patch);
 795
 796     if (patch >= 0)
 797     {
 798         // Modify patch information for this boundary edge
 799         edgePatchSizes_[patch]++;
 800
 801         for (label i = (patch + 1); i < nPatches_; i++)
 802         {
 803             edgePatchStarts_[i]++;
 804         }
 805     }
 806     else
 807     {
 808         // Increment the number of internal edges, and subsequent patch-starts
 809         nInternalEdges_++;
 810
 811         for (label i = 0; i < nPatches_; i++)
 812         {
 813             edgePatchStarts_[i]++;
 814         }
 815     }
 816
 817     // Size-up the pointEdges list
 818     if (!twoDMesh_)
 819     {
 820         meshOps::sizeUpList(newEdgeIndex, pointEdges_[newEdge[0]]);
 821         meshOps::sizeUpList(newEdgeIndex, pointEdges_[newEdge[1]]);
 822     }
 823
 824     // Increment the total edge count
 825     nEdges_++;
 826
 827     return newEdgeIndex;
 828 }
 829
 830
 831 // Remove the specified edge from the mesh
 832 void dynamicTopoFvMesh::removeEdge
 833 (
 834     const label eIndex
 835 )
 836 {
 837     if (!twoDMesh_)
 838     {
 839         const edge& rEdge = edges_[eIndex];
 840
 841         // Size-down the pointEdges list
 842         if (pointEdges_[rEdge[0]].size())
 843         {
 844             meshOps::sizeDownList(eIndex, pointEdges_[rEdge[0]]);
 845         }
 846
 847         if (pointEdges_[rEdge[1]].size())
 848         {
 849             meshOps::sizeDownList(eIndex, pointEdges_[rEdge[1]]);
 850         }
 851
 852         // Entity won't be removed from the stack for efficiency
 853         // It will be discarded on access instead.
 854
 855         // Update coupled edge maps, if necessary.
 856         forAll(patchCoupling_, patchI)
 857         {
 858             if (patchCoupling_(patchI))
 859             {
 860                 const coupleMap& cMap = patchCoupling_[patchI].map();
 861
 862                 cMap.removeMaster(coupleMap::EDGE, eIndex);
 863                 cMap.removeSlave(coupleMap::EDGE, eIndex);
 864             }
 865         }
 866     }
 867
 868     // Identify the patch for this edge
 869     label patch = whichEdgePatch(eIndex);
 870
 871     if (debug > 2)
 872     {
 873         Pout<< "Removing edge: "
 874             << eIndex << ": "
 875             << edges_[eIndex];
 876
 877         const polyBoundaryMesh& boundary = boundaryMesh();
 878
 879         Pout<< " Patch: ";
 880
 881         if (patch == -1)
 882         {
 883             Pout<< "Internal" << endl;
 884         }
 885         else
 886         if (patch < boundary.size())
 887         {
 888             Pout<< boundary[patch].name() << endl;
 889         }
 890         else
 891         {
 892             Pout<< " New patch: " << patch << endl;
 893         }
 894     }
 895
 896     // Invalidate edge and edgeFaces
 897     edges_[eIndex] = edge(-1, -1);
 898     edgeFaces_[eIndex].clear();
 899
 900     if (patch >= 0)
 901     {
 902         // Modify patch information for this boundary edge
 903         edgePatchSizes_[patch]--;
 904
 905         for (label i = (patch + 1); i < nPatches_; i++)
 906         {
 907             edgePatchStarts_[i]--;
 908         }
 909     }
 910     else
 911     {
 912         // Decrement the internal edge count, and subsequent patch-starts
 913         nInternalEdges_--;
 914
 915         forAll(edgePatchStarts_, patchI)
 916         {
 917             edgePatchStarts_[patchI]--;
 918         }
 919     }
 920
 921     // Update reverse edge-map, but only if this is an edge that existed
 922     // at time [n]. Added edges which are deleted during the topology change
 923     // needn't be updated.
 924     if (eIndex < nOldEdges_)
 925     {
 926         reverseEdgeMap_[eIndex] = -1;
 927     }
 928     else
 929     {
 930         // Store this information for the reOrdering stage
 931         deletedEdges_.insert(eIndex);
 932     }
 933
 934     // Check and remove from the list of added edge patches
 935     Map<label>::iterator it = addedEdgePatches_.find(eIndex);
 936
 937     if (it != addedEdgePatches_.end())
 938     {
 939         addedEdgePatches_.erase(it);
 940     }
 941
 942     // Decrement the total edge-count
 943     nEdges_--;
 944 }
 945
 946
 947 // Insert the specified point to the mesh
 948 label dynamicTopoFvMesh::insertPoint
 949 (
 950     const point& newPoint,
 951     const point& oldPoint,
 952     const labelList& mapPoints,
 953     const label zoneID
 954 )
 955 {
 956     // Add a new point to the end of the list
 957     label newPointIndex = points_.size();
 958
 959     points_.append(newPoint);
 960     oldPoints_.append(oldPoint);
 961
 962     if (debug > 2)
 963     {
 964         Pout<< "Inserting new point: "
 965             << newPointIndex << ": "
 966             << newPoint
 967             << " and old point: "
 968             << oldPoint
 969             << "  Mapped from: "
 970             << mapPoints << endl;
 971     }
 972
 973     // Make a pointsFromPoints entry
 974     meshOps::sizeUpList
 975     (
 976         objectMap(newPointIndex, mapPoints),
 977         pointsFromPoints_
 978     );
 979
 980     // Add an empty entry to pointEdges as well.
 981     // This entry can be sized-up appropriately at a later stage.
 982     if (!twoDMesh_)
 983     {
 984         pointEdges_.append(labelList(0));
 985     }
 986
 987     // Add to the zone if necessary
 988     if (zoneID >= 0)
 989     {
 990         addedPointZones_.insert(newPointIndex, zoneID);
 991     }
 992
 993     nPoints_++;
 994
 995     return newPointIndex;
 996 }
 997
 998
 999 // Remove the specified point from the mesh
1000 void dynamicTopoFvMesh::removePoint
1001 (
1002     const label pIndex
1003 )
1004 {
1005     if (debug > 2)
1006     {
1007         Pout<< "Removing point: " << pIndex
1008             << " :: new: " << points_[pIndex]
1009             << " :: old: " << oldPoints_[pIndex]
1010             << endl;
1011     }
1012
1013     // Remove the point
1014     // (or just make sure that it's never used anywhere else)
1015     // points_[pIndex] = point();
1016     // oldPoints_[pIndex] = point();
1017
1018     // Remove pointEdges as well
1019     if (!twoDMesh_)
1020     {
1021         pointEdges_[pIndex].clear();
1022     }
1023
1024     // Update the reverse point map
1025     if (pIndex < nOldPoints_)
1026     {
1027         reversePointMap_[pIndex] = -1;
1028     }
1029     else
1030     {
1031         deletedPoints_.insert(pIndex);
1032     }
1033
1034     // Check if this point was added to a zone
1035     Map<label>::iterator pzit = addedPointZones_.find(pIndex);
1036
1037     if (pzit != addedPointZones_.end())
1038     {
1039         addedPointZones_.erase(pzit);
1040     }
1041
1042     // Update coupled point maps, if necessary.
1043     forAll(patchCoupling_, patchI)
1044     {
1045         if (patchCoupling_(patchI))
1046         {
1047             // Obtain references
1048             const coupleMap & cMap = patchCoupling_[patchI].map();
1049
1050             Map<label>& pointMap = cMap.entityMap(coupleMap::POINT);
1051             Map<label>& rPointMap = cMap.reverseEntityMap(coupleMap::POINT);
1052
1053             Map<label>::iterator pmit = pointMap.find(pIndex);
1054
1055             if (pmit != pointMap.end())
1056             {
1057                 // Erase the reverse map first
1058                 rPointMap.erase(pmit());
1059
1060                 // Update pointMap
1061                 pointMap.erase(pmit);
1062             }
1063         }
1064     }
1065
1066     // Check if the point was added in the current morph, and delete
1067     forAll(pointsFromPoints_, indexI)
1068     {
1069         if (pointsFromPoints_[indexI].index() == pIndex)
1070         {
1071             // Remove entry from the list
1072             meshOps::removeIndex(indexI, pointsFromPoints_);
1073             break;
1074         }
1075     }
1076
1077     // Decrement the total point-count
1078     nPoints_--;
1079 }
1080
1081
1082 // Utility method to build vertexHull for an edge [3D].
1083 // Assumes that edgeFaces information is consistent.
1084 void dynamicTopoFvMesh::buildVertexHull
1085 (
1086     const label eIndex,
1087     labelList& vertexHull,
1088     const label checkIndex
1089 ) const
1090 {
1091     bool found = false;
1092     label faceIndex = -1, cellIndex = -1;
1093     label otherPoint = -1, nextPoint = -1;
1094     label failMode = 0;
1095
1096     // Obtain references
1097     const edge& edgeToCheck = edges_[eIndex];
1098     const labelList& eFaces = edgeFaces_[eIndex];
1099
1100     // Re-size the list first
1101     vertexHull.clear();
1102     vertexHull.setSize(eFaces.size(), -1);
1103
1104     if (whichEdgePatch(eIndex) == -1)
1105     {
1106         // Internal edge.
1107         // Pick the first face and start with that
1108         faceIndex = eFaces[0];
1109     }
1110     else
1111     {
1112         // Need to find a properly oriented start-face
1113         forAll(eFaces, faceI)
1114         {
1115             if (whichPatch(eFaces[faceI]) > -1)
1116             {
1117                 meshOps::findIsolatedPoint
1118                 (
1119                     faces_[eFaces[faceI]],
1120                     edgeToCheck,
1121                     otherPoint,
1122                     nextPoint
1123                 );
1124
1125                 if (nextPoint == edgeToCheck[checkIndex])
1126                 {
1127                     faceIndex = eFaces[faceI];
1128                     break;
1129                 }
1130             }
1131         }
1132     }
1133
1134     if (faceIndex == -1)
1135     {
1136         // Define failure-mode
1137         failMode = 1;
1138     }
1139     else
1140     {
1141         // Shuffle vertices to appear in CCW order
1142         forAll(vertexHull, indexI)
1143         {
1144             meshOps::findIsolatedPoint
1145             (
1146                 faces_[faceIndex],
1147                 edgeToCheck,
1148                 otherPoint,
1149                 nextPoint
1150             );
1151
1152             // Add the isolated point
1153             vertexHull[indexI] = otherPoint;
1154
1155             // Figure out how this edge is oriented.
1156             if (nextPoint == edgeToCheck[checkIndex])
1157             {
1158                 // Counter-clockwise. Pick the owner.
1159                 cellIndex = owner_[faceIndex];
1160             }
1161             else
1162             if (whichPatch(faceIndex) == -1)
1163             {
1164                 // Clockwise. Pick the neighbour.
1165                 cellIndex = neighbour_[faceIndex];
1166             }
1167             else
1168             if (indexI != (vertexHull.size() - 1))
1169             {
1170                 // This could be a pinched manifold edge
1171                 // (situation with more than two boundary faces)
1172                 // Pick another (properly oriented) boundary face.
1173                 faceIndex = -1;
1174
1175                 forAll(eFaces, faceI)
1176                 {
1177                     if (whichPatch(eFaces[faceI]) > -1)
1178                     {
1179                         meshOps::findIsolatedPoint
1180                         (
1181                             faces_[eFaces[faceI]],
1182                             edgeToCheck,
1183                             otherPoint,
1184                             nextPoint
1185                         );
1186
1187                         if
1188                         (
1189                             (nextPoint == edgeToCheck[checkIndex]) &&
1190                             (findIndex(vertexHull, otherPoint) == -1)
1191                         )
1192                         {
1193                             faceIndex = eFaces[faceI];
1194                             break;
1195                         }
1196                     }
1197                 }
1198
1199                 if (faceIndex == -1)
1200                 {
1201                     failMode = 2;
1202                     break;
1203                 }
1204                 else
1205                 {
1206                     continue;
1207                 }
1208             }
1209             else
1210             {
1211                 // Looks like we've hit the last boundary face. Break out.
1212                 break;
1213             }
1214
1215             const cell& cellToCheck = cells_[cellIndex];
1216
1217             found = false;
1218
1219             // Assuming tet-cells,
1220             // Loop through edgeFaces and get the next face
1221             forAll(eFaces, faceI)
1222             {
1223                 if
1224                 (
1225                     eFaces[faceI] != faceIndex
1226                  && eFaces[faceI] == cellToCheck[0]
1227                 )
1228                 {
1229                     faceIndex = cellToCheck[0];
1230                     found = true; break;
1231                 }
1232
1233                 if
1234                 (
1235                     eFaces[faceI] != faceIndex
1236                  && eFaces[faceI] == cellToCheck[1]
1237                 )
1238                 {
1239                     faceIndex = cellToCheck[1];
1240                     found = true; break;
1241                 }
1242
1243                 if
1244                 (
1245                     eFaces[faceI] != faceIndex
1246                  && eFaces[faceI] == cellToCheck[2]
1247                 )
1248                 {
1249                     faceIndex = cellToCheck[2];
1250                     found = true; break;
1251                 }
1252
1253                 if
1254                 (
1255                     eFaces[faceI] != faceIndex
1256                  && eFaces[faceI] == cellToCheck[3]
1257                 )
1258                 {
1259                     faceIndex = cellToCheck[3];
1260                     found = true; break;
1261                 }
1262             }
1263
1264             if (!found)
1265             {
1266                 failMode = 3;
1267                 break;
1268             }
1269         }
1270     }
1271
1272     if (debug > 2 && !failMode)
1273     {
1274         // Check for invalid indices
1275         if (findIndex(vertexHull, -1) > -1)
1276         {
1277             failMode = 4;
1278         }
1279
1280         // Check for duplicate labels
1281         if (!failMode)
1282         {
1283             labelHashSet uniquePoints;
1284
1285             forAll(vertexHull, pointI)
1286             {
1287                 bool inserted = uniquePoints.insert(vertexHull[pointI]);
1288
1289                 if (!inserted)
1290                 {
1291                     Pout<< " vertexHull for edge: "
1292                         << eIndex << "::" << edgeToCheck
1293                         << " contains identical vertex labels: "
1294                         << vertexHull << endl;
1295
1296                     failMode = 5;
1297                 }
1298             }
1299         }
1300     }
1301
1302     if (failMode)
1303     {
1304         // Prepare edgeCells
1305         DynamicList<label> eCells(10);
1306
1307         forAll(eFaces, faceI)
1308         {
1309             label own = owner_[eFaces[faceI]];
1310             label nei = neighbour_[eFaces[faceI]];
1311
1312             if (findIndex(eCells, own) == -1)
1313             {
1314                 eCells.append(own);
1315             }
1316
1317             if (nei > -1)
1318             {
1319                 if (findIndex(eCells, nei) == -1)
1320                 {
1321                     eCells.append(nei);
1322                 }
1323             }
1324         }
1325
1326         // Write out for post-processing
1327         label eI = eIndex, epI = whichEdgePatch(eIndex);
1328         word epName(epI < 0 ? "Internal" : boundaryMesh()[epI].name());
1329
1330         writeVTK("vRingEdge_" + Foam::name(eI), eIndex, 1, false, true);
1331         writeVTK("vRingEdgeFaces_" + Foam::name(eI), eFaces, 2, false, true);
1332         writeVTK("vRingEdgeCells_" + Foam::name(eI), eCells, 3, false, true);
1333
1334         Pout<< "edgeFaces: " << endl;
1335
1336         forAll(eFaces, faceI)
1337         {
1338             label fpI = whichPatch(eFaces[faceI]);
1339             word fpName(fpI < 0 ? "Internal" : boundaryMesh()[fpI].name());
1340
1341             Pout<< " Face: " << eFaces[faceI]
1342                 << ":: " << faces_[eFaces[faceI]]
1343                 << " Owner: " << owner_[eFaces[faceI]]
1344                 << " Neighbour: " << neighbour_[eFaces[faceI]]
1345                 << " Patch: " << fpName
1346                 << endl;
1347         }
1348
1349         FatalErrorIn
1350         (
1351             "\n"
1352             "void dynamicTopoFvMesh::buildVertexHull\n"
1353             "(\n"
1354             "    const label eIndex,\n"
1355             "    labelList& vertexHull,\n"
1356             "    const label checkIndex\n"
1357             ")\n"
1358         )
1359             << " Failed to determine a vertex ring. " << nl
1360             << " Failure mode: " << failMode << nl
1361             << " Edge: " << eIndex << ":: " << edgeToCheck << nl
1362             << " edgeFaces: " << eFaces << nl
1363             << " Patch: " << epName << nl
1364             << " Current vertexHull: " << vertexHull
1365             << abort(FatalError);
1366     }
1367 }
1368
1369
1370 void dynamicTopoFvMesh::handleMeshSlicing()
1371 {
1372     if (slicePairs_.empty())
1373     {
1374         return;
1375     }
1376
1377     if (lengthEstimator().holdOff())
1378     {
1379         // Clear out data.
1380         slicePairs_.clear();
1381         lengthEstimator().clearBoxes();
1382
1383         // Hold-off mesh slicing for a few time-steps.
1384         lengthEstimator().decrementHoldOff();
1385
1386         return;
1387     }
1388
1389     Info<< "Slicing Mesh...";
1390
1391     // Loop through candidates and weed-out invalid points
1392     forAll(slicePairs_, pairI)
1393     {
1394         const labelPair& pairToCheck = slicePairs_[pairI];
1395
1396         bool available = true;
1397
1398         if (twoDMesh_)
1399         {
1400             forAll(pairToCheck, indexI)
1401             {
1402                 if (deletedFaces_.found(pairToCheck[indexI]))
1403                 {
1404                     available = false;
1405                     break;
1406                 }
1407
1408                 if (pairToCheck[indexI] < nOldFaces_)
1409                 {
1410                     if (reverseFaceMap_[pairToCheck[indexI]] == -1)
1411                     {
1412                         available = false;
1413                         break;
1414                     }
1415                 }
1416             }
1417         }
1418         else
1419         {
1420             forAll(pairToCheck, indexI)
1421             {
1422                 if (deletedPoints_.found(pairToCheck[indexI]))
1423                 {
1424                     available = false;
1425                     break;
1426                 }
1427
1428                 if (pairToCheck[indexI] < nOldPoints_)
1429                 {
1430                     if (reversePointMap_[pairToCheck[indexI]] == -1)
1431                     {
1432                         available = false;
1433                         break;
1434                     }
1435                 }
1436             }
1437         }
1438
1439         if (available)
1440         {
1441             // Slice the mesh at this point.
1442             sliceMesh(pairToCheck);
1443         }
1444     }
1445
1446     Info<< "Done." << endl;
1447
1448     // Clear out data.
1449     slicePairs_.clear();
1450     lengthEstimator().clearBoxes();
1451
1452     // Set the sliceHoldOff value
1453     lengthEstimator().setHoldOff(50);
1454 }
1455
1456
1457 // Handle layer addition / removal events
1458 void dynamicTopoFvMesh::handleLayerAdditionRemoval()
1459 {
1460     const dictionary& meshSubDict = dict_.subDict("dynamicTopoFvMesh");
1461
1462     if (meshSubDict.found("layering") || mandatory_)
1463     {
1464         // Bail out if no entries are present
1465         if (meshSubDict.subDict("layering").empty())
1466         {
1467             return;
1468         }
1469     }
1470     else
1471     {
1472         // No dictionary present.
1473         return;
1474     }
1475
1476     // Fetch layering patches
1477     const polyBoundaryMesh& boundary = boundaryMesh();
1478     const dictionary& layeringDict = meshSubDict.subDict("layering");
1479
1480     wordList layeringPatches = layeringDict.toc();
1481
1482     forAll(layeringPatches, wordI)
1483     {
1484         word pName = layeringPatches[wordI];
1485         label patchID = boundary.findPatchID(pName);
1486
1487         // If this patch has no faces, move on
1488         if (boundary[patchID].empty())
1489         {
1490             continue;
1491         }
1492
1493         // Use first face to determine layer thickness
1494         scalar magSf = mag(boundary[patchID].faceAreas()[0]);
1495         scalar Vc = cellVolumes()[boundary[patchID].faceCells()[0]];
1496
1497         // Define thickness
1498         scalar thickness = (Vc / (magSf + VSMALL));
1499
1500         // Fetch min / max thickness
1501         scalar minThickness =
1502         (
1503             readScalar(layeringDict.subDict(pName).lookup("minThickness"))
1504         );
1505
1506         scalar maxThickness =
1507         (
1508             readScalar(layeringDict.subDict(pName).lookup("maxThickness"))
1509         );
1510
1511         if (debug)
1512         {
1513             Pout<< " Patch: " << pName << nl
1514                 << " Face area: " << magSf << nl
1515                 << " Cell volume: " << Vc << nl
1516                 << " Layer thickness: " << thickness << nl
1517                 << " Min thickness: " << minThickness << nl
1518                 << " Max thickness: " << maxThickness << nl
1519                 << endl;
1520         }
1521
1522         if (thickness > maxThickness)
1523         {
1524             // Add cell layer above patch
1525             addCellLayer(patchID);
1526         }
1527         else
1528         if (thickness < minThickness)
1529         {
1530             // Remove cell layer above patch
1531             removeCellLayer(patchID);
1532         }
1533     }
1534 }
1535
1536
1537 // Test an edge / face for proximity with other faces on proximity patches
1538 // and return the scalar distance to an oppositely-oriented face.
1539 scalar dynamicTopoFvMesh::testProximity
1540 (
1541     const label index,
1542     label& proximityFace
1543 ) const
1544 {
1545     scalar proxDistance = GREAT, testStep = 0.0;
1546     vector gCentre = vector::zero, gNormal = vector::zero;
1547
1548     if (twoDMesh_)
1549     {
1550         // Obtain the face-normal.
1551         gNormal = faces_[index].normal(points_);
1552
1553         // Obtain the face centre.
1554         gCentre = faces_[index].centre(points_);
1555
1556         // Fetch the edge
1557         const edge& edgeToCheck = edges_[getTriBoundaryEdge(index)];
1558
1559         // Calculate a test step-size
1560         testStep =
1561         (
1562             linePointRef
1563             (
1564                 points_[edgeToCheck.start()],
1565                 points_[edgeToCheck.end()]
1566             ).mag()
1567         );
1568     }
1569     else
1570     {
1571         const edge& edgeToCheck = edges_[index];
1572         const labelList& eFaces = edgeFaces_[index];
1573
1574         linePointRef lpr
1575         (
1576             points_[edgeToCheck.start()],
1577             points_[edgeToCheck.end()]
1578         );
1579
1580         // Obtain the edge centre.
1581         gCentre = lpr.centre();
1582
1583         // Obtain the edge-normal
1584         forAll(eFaces, faceI)
1585         {
1586             if (neighbour_[eFaces[faceI]] == -1)
1587             {
1588                 // Obtain the normal.
1589                 gNormal += faces_[eFaces[faceI]].normal(points_);
1590             }
1591         }
1592
1593         // Calculate a test step-size
1594         testStep = lpr.mag();
1595     }
1596
1597     // Normalize
1598     gNormal /= (mag(gNormal) + VSMALL);
1599
1600     // Test for proximity, and mark slice-pairs
1601     // is the distance is below threshold.
1602     if
1603     (
1604         lengthEstimator().testProximity
1605         (
1606             gCentre,
1607             gNormal,
1608             testStep,
1609             proximityFace,
1610             proxDistance
1611         )
1612     )
1613     {
1614         labelPair proxPoints(-1, -1);
1615         bool foundPoint = false;
1616
1617         if (twoDMesh_)
1618         {
1619             proxPoints.first() = index;
1620             proxPoints.second() = proximityFace;
1621
1622             if
1623             (
1624                 (faces_[index].size() == 4) &&
1625                 (polyMesh::faces()[proximityFace].size() == 4)
1626             )
1627             {
1628                 foundPoint = true;
1629             }
1630         }
1631         else
1632         {
1633             const edge& thisEdge = edges_[index];
1634             const face& proxFace = polyMesh::faces()[proximityFace];
1635
1636             // Check if any points on this face are still around.
1637             // If yes, mark one of them as the end point
1638             // for Dijkstra's algorithm. The start point will be a point
1639             // on this edge.
1640             proxPoints.first() = thisEdge[0];
1641
1642             forAll(proxFace, pointI)
1643             {
1644                 if (reversePointMap_[proxFace[pointI]] != -1)
1645                 {
1646                     proxPoints.second() = proxFace[pointI];
1647                     foundPoint = true;
1648                     break;
1649                 }
1650             }
1651         }
1652
1653         if (foundPoint)
1654         {
1655             // Lock the edge mutex
1656             entityMutex_[1].lock();
1657
1658             label newSize = slicePairs_.size() + 1;
1659
1660             // Const-cast slicePairs for resize
1661             List<labelPair>& sP = const_cast<List<labelPair>&>(slicePairs_);
1662
1663             // Add this entry as a candidate for mesh slicing.
1664             sP.setSize(newSize, proxPoints);
1665
1666             // Unlock the edge mutex
1667             entityMutex_[1].unlock();
1668         }
1669     }
1670
1671     return proxDistance;
1672 }
1673
1674
1675 // Calculate the edge length-scale for the mesh
1676 void dynamicTopoFvMesh::calculateLengthScale(bool dump)
1677 {
1678     if (!edgeRefinement_)
1679     {
1680         return;
1681     }
1682
1683     Switch dumpLengthScale(false);
1684
1685     const dictionary& meshDict = dict_.subDict("dynamicTopoFvMesh");
1686
1687     if (meshDict.found("dumpLengthScale") || mandatory_)
1688     {
1689         dumpLengthScale = readBool(meshDict.lookup("dumpLengthScale"));
1690     }
1691
1692     autoPtr<volScalarField> lsfPtr(NULL);
1693
1694     if (dumpLengthScale && time().outputTime() && dump)
1695     {
1696         lsfPtr.set
1697         (
1698             new volScalarField
1699             (
1700                 IOobject
1701                 (
1702                     "lengthScale",
1703                     time().timeName(),
1704                     *this,
1705                     IOobject::NO_READ,
1706                     IOobject::NO_WRITE,
1707                     false
1708                 ),
1709                 *this,
1710                 dimensionedScalar("scalar", dimLength, 0)
1711             )
1712         );
1713     }
1714
1715     // Bail-out if a dumping was not requested in dictionary.
1716     if (dump && !dumpLengthScale)
1717     {
1718         return;
1719     }
1720
1721     // Size the field and calculate length-scale
1722     lengthScale_.setSize(nCells_, 0.0);
1723
1724     lengthEstimator().calculateLengthScale(lengthScale_);
1725
1726     // Check if length-scale is to be dumped to disk.
1727     if (dumpLengthScale && time().outputTime() && dump)
1728     {
1729         // Obtain length-scale values from the mesh
1730         lsfPtr->internalField() = lengthScale_;
1731
1732         lsfPtr->write();
1733     }
1734 }
1735
1736
1737 // Read optional dictionary parameters
1738 void dynamicTopoFvMesh::readOptionalParameters(bool reRead)
1739 {
1740     // Read from disk
1741     dict_.readIfModified();
1742
1743     const dictionary& meshSubDict = dict_.subDict("dynamicTopoFvMesh");
1744
1745     // Enable/disable run-time debug level
1746     if (meshSubDict.found("debug") || mandatory_)
1747     {
1748         // Look for a processor-specific debug flag
1749         if (Pstream::parRun() && meshSubDict.found("parDebugProcs"))
1750         {
1751             labelList procs(meshSubDict.lookup("parDebugProcs"));
1752
1753             forAll(procs, procI)
1754             {
1755                 if (Pstream::myProcNo() == procs[procI])
1756                 {
1757                     debug = readLabel(meshSubDict.lookup("debug"));
1758                 }
1759             }
1760         }
1761         else
1762         {
1763             debug = readLabel(meshSubDict.lookup("debug"));
1764         }
1765     }
1766     else
1767     {
1768         debug = 0;
1769     }
1770
1771     // Set debug option for underlying classes as well.
1772     if (debug > 2)
1773     {
1774         fvMesh::debug = true;
1775         polyMesh::debug = true;
1776     }
1777
1778     // Re-read edge-refinement options, if necessary
1779     if (edgeRefinement_ && reRead)
1780     {
1781         lengthEstimator().readRefinementOptions
1782         (
1783             meshSubDict.subDict("refinementOptions"),
1784             true,
1785             mandatory_
1786         );
1787     }
1788
1789     // Read re-mesh interval
1790     if (meshSubDict.found("interval") || mandatory_)
1791     {
1792         interval_ = readLabel(meshSubDict.lookup("interval"));
1793     }
1794     else
1795     {
1796         interval_ = 1;
1797     }
1798
1799     // Check if an external mesh-motion solver is used
1800     if (meshSubDict.found("loadMotionSolver") || mandatory_)
1801     {
1802         loadMotionSolver_.readIfPresent("loadMotionSolver", meshSubDict);
1803     }
1804
1805     // Update bandwidth reduction switch
1806     if (meshSubDict.found("bandwidthReduction") || mandatory_)
1807     {
1808         bandWidthReduction_.readIfPresent("bandwidthReduction", meshSubDict);
1809     }
1810
1811     // Update threshold for sliver cells
1812     if (meshSubDict.found("sliverThreshold") || mandatory_)
1813     {
1814         sliverThreshold_ = readScalar(meshSubDict.lookup("sliverThreshold"));
1815
1816         if (sliverThreshold_ > 1.0 || sliverThreshold_ < 0.0)
1817         {
1818             FatalErrorIn("void dynamicTopoFvMesh::readOptionalParameters()")
1819                 << " Sliver threshold out of range [0..1]"
1820                 << abort(FatalError);
1821         }
1822     }
1823
1824     // Update limit for max number of bisections / collapses
1825     if (meshSubDict.found("maxModifications") || mandatory_)
1826     {
1827         maxModifications_ = readLabel(meshSubDict.lookup("maxModifications"));
1828     }
1829
1830     // Update limit for swap on curved surfaces
1831     if (meshSubDict.found("swapDeviation") || mandatory_)
1832     {
1833         swapDeviation_ = readScalar(meshSubDict.lookup("swapDeviation"));
1834
1835         if (swapDeviation_ > 1.0 || swapDeviation_ < 0.0)
1836         {
1837             FatalErrorIn("void dynamicTopoFvMesh::readOptionalParameters()")
1838                 << " Swap deviation out of range [0..1]"
1839                 << abort(FatalError);
1840         }
1841     }
1842
1843     // For tetrahedral meshes...
1844     if (!twoDMesh_)
1845     {
1846         // Check if swapping is to be avoided on any patches
1847         if (meshSubDict.found("noSwapPatches") || mandatory_)
1848         {
1849             wordList noSwapPatches =
1850             (
1851                 meshSubDict.subDict("noSwapPatches").toc()
1852             );
1853
1854             noSwapPatchIDs_.setSize(noSwapPatches.size());
1855
1856             label indexI = 0;
1857
1858             forAll(noSwapPatches, wordI)
1859             {
1860                 const word& patchName = noSwapPatches[wordI];
1861
1862                 noSwapPatchIDs_[indexI++] =
1863                 (
1864                     boundaryMesh().findPatchID(patchName)
1865                 );
1866             }
1867         }
1868
1869         // Check if a limit has been imposed on maxTetsPerEdge
1870         if (meshSubDict.found("maxTetsPerEdge") || mandatory_)
1871         {
1872             maxTetsPerEdge_ = readLabel(meshSubDict.lookup("maxTetsPerEdge"));
1873         }
1874         else
1875         {
1876             maxTetsPerEdge_ = 7;
1877         }
1878
1879         // Check if programming tables can be resized at runtime
1880         if (meshSubDict.found("allowTableResize") || mandatory_)
1881         {
1882             allowTableResize_ =
1883             (
1884                 readBool(meshSubDict.lookup("allowTableResize"))
1885             );
1886         }
1887         else
1888         {
1889             allowTableResize_ = false;
1890         }
1891     }
1892
1893     // Check for load-balancing in parallel
1894     if (reRead && (meshSubDict.found("loadBalancing") || mandatory_))
1895     {
1896         // Execute balancing
1897         executeLoadBalancing(meshSubDict.subDict("loadBalancing"));
1898     }
1899 }
1900
1901
1902 // Initialize edge related connectivity lists
1903 void dynamicTopoFvMesh::initEdges()
1904 {
1905     // Initialize eMesh, and copy to local lists
1906     eMeshPtr_.set(new eMesh(*this));
1907
1908     // Obtain information
1909     nEdges_ = eMeshPtr_->nEdges();
1910     nInternalEdges_ = eMeshPtr_->nInternalEdges();
1911     edgePatchSizes_ = eMeshPtr_->boundary().patchSizes();
1912     edgePatchStarts_ = eMeshPtr_->boundary().patchStarts();
1913
1914     // Set old edge information
1915     nOldEdges_ = nEdges_;
1916     nOldInternalEdges_ = nInternalEdges_;
1917     oldEdgePatchSizes_ = edgePatchSizes_;
1918     oldEdgePatchStarts_ = edgePatchStarts_;
1919
1920     // Set local lists with edge connectivity information
1921     edges_ = eMeshPtr_->edges();
1922     edgeFaces_ = eMeshPtr_->edgeFaces();
1923     faceEdges_ = eMeshPtr_->faceEdges();
1924
1925     if (!twoDMesh_)
1926     {
1927         // Invert edges to obtain pointEdges
1928         pointEdges_ = invertManyToMany<edge, labelList>(nPoints_, edges_);
1929     }
1930
1931     // Clear out unwanted eMesh connectivity
1932     eMeshPtr_->clearOut();
1933
1934     // Read coupled patch information from dictionary.
1935     readCoupledPatches();
1936 }
1937
1938
1939 // Load the mesh-quality metric from the library
1940 void dynamicTopoFvMesh::loadMetric()
1941 {
1942     if (twoDMesh_)
1943     {
1944         return;
1945     }
1946
1947     // Specify the dictionary we would be looking in...
1948     const dictionary& meshDict = dict_.subDict("dynamicTopoFvMesh");
1949
1950     // Select an appropriate metric
1951     tetMetric_ = tetMetric::New(meshDict, meshDict.lookup("tetMetric"));
1952 }
1953
1954
1955 // Load the mesh-motion solver
1956 void dynamicTopoFvMesh::loadMotionSolver()
1957 {
1958     if (motionSolver_.valid())
1959     {
1960         FatalErrorIn
1961         (
1962             "void dynamicTopoFvMesh::loadMotionSolver() "
1963         ) << nl << " Motion solver already loaded. "
1964           << abort(FatalError);
1965     }
1966     else
1967     if (dict_.found("solver") && loadMotionSolver_)
1968     {
1969         motionSolver_ = motionSolver::New(*this);
1970     }
1971 }
1972
1973
1974 // Load the field mapper
1975 void dynamicTopoFvMesh::loadFieldMapper()
1976 {
1977     if (mapper_.valid())
1978     {
1979         FatalErrorIn
1980         (
1981             "void dynamicTopoFvMesh::loadFieldMapper() "
1982         ) << nl << " Field mapper already loaded. "
1983           << abort(FatalError);
1984     }
1985     else
1986     {
1987         mapper_.set(new topoMapper(*this, dict_));
1988     }
1989 }
1990
1991
1992 // Load the length scale estimator
1993 void dynamicTopoFvMesh::loadLengthScaleEstimator()
1994 {
1995     if (edgeRefinement_)
1996     {
1997         if (lengthEstimator_.valid())
1998         {
1999             FatalErrorIn
2000             (
2001                 "void dynamicTopoFvMesh::loadLengthScaleEstimator() "
2002             ) << nl << " Length scale estimator already loaded. "
2003               << abort(FatalError);
2004         }
2005         else
2006         {
2007             lengthEstimator_.set
2008             (
2009                 new lengthScaleEstimator(*this)
2010             );
2011         }
2012
2013         // Read options
2014         lengthEstimator().readRefinementOptions
2015         (
2016             dict_.subDict("dynamicTopoFvMesh").subDict("refinementOptions"),
2017             false,
2018             mandatory_
2019         );
2020
2021         // Set coupled patch options, if available
2022         if (dict_.found("coupledPatches") || mandatory_)
2023         {
2024             lengthEstimator().setCoupledPatches
2025             (
2026                 dict_.subDict("coupledPatches")
2027             );
2028         }
2029     }
2030 }
2031
2032
2033 // Initialize the threading environment.
2034 //  - Provides an override option to avoid reading from the dictionary.
2035 void dynamicTopoFvMesh::initializeThreadingEnvironment
2036 (
2037     const label specThreads
2038 )
2039 {
2040     // Initialize an IOobject for the IOmultiThreader object
2041     IOobject io
2042     (
2043         "threader",
2044         this->time().timeName(),
2045         (*this),
2046         IOobject::NO_READ,
2047         IOobject::NO_WRITE,
2048         true
2049     );
2050
2051     if (specThreads > 0)
2052     {
2053         threader_.set(new IOmultiThreader(io, specThreads));
2054     }
2055     else
2056     {
2057         if (dict_.subDict("dynamicTopoFvMesh").found("threads") || mandatory_)
2058         {
2059             threader_.set
2060             (
2061                 new IOmultiThreader
2062                 (
2063                     io,
2064                     readLabel
2065                     (
2066                         dict_.subDict("dynamicTopoFvMesh").lookup("threads")
2067                     )
2068                 )
2069             );
2070         }
2071         else
2072         {
2073             threader_.set(new IOmultiThreader(io, 1));
2074         }
2075     }
2076
2077     // Get the number of threads and allocate threadHandlers
2078     label nThreads = threader_->getNumThreads();
2079
2080     if (nThreads == 1)
2081     {
2082         handlerPtr_.setSize(1);
2083
2084         handlerPtr_.set
2085         (
2086             0,
2087             new meshHandler(*this, threader())
2088         );
2089
2090         handlerPtr_[0].setMaster();
2091
2092         // Size the stacks
2093         entityStack_.setSize(1);
2094     }
2095     else
2096     {
2097         // Index '0' is master, rest are slaves
2098         handlerPtr_.setSize(nThreads + 1);
2099
2100         // Size the stacks
2101         entityStack_.setSize(nThreads + 1);
2102
2103         forAll(handlerPtr_, threadI)
2104         {
2105             handlerPtr_.set
2106             (
2107                 threadI,
2108                 new meshHandler(*this, threader())
2109             );
2110
2111             if (threadI == 0)
2112             {
2113                 handlerPtr_[0].setMaster();
2114             }
2115             else
2116             {
2117                 handlerPtr_[threadI].setID(threader_->getID(threadI - 1));
2118                 handlerPtr_[threadI].setSlave();
2119             }
2120         }
2121     }
2122 }
2123
2124
2125 // 2D Edge-swapping engine
2126 void dynamicTopoFvMesh::swap2DEdges(void *argument)
2127 {
2128     // Recast the argument
2129     meshHandler *thread = static_cast<meshHandler*>(argument);
2130
2131     if (thread->slave())
2132     {
2133         thread->sendSignal(meshHandler::START);
2134     }
2135
2136     dynamicTopoFvMesh& mesh = thread->reference();
2137
2138     // Figure out which thread this is...
2139     label tIndex = mesh.self();
2140
2141     // Set the timer
2142     clockTime sTimer;
2143
2144     bool reported = false;
2145     label stackSize = mesh.stack(tIndex).size();
2146     scalar interval = mesh.reportInterval(), oIndex = 0.0, nIndex = 0.0;
2147
2148     oIndex = ::floor(sTimer.elapsedTime() / interval);
2149
2150     // Pick items off the stack
2151     while (!mesh.stack(tIndex).empty())
2152     {
2153         // Report progress
2154         if (thread->master())
2155         {
2156             // Update the index, if its changed
2157             nIndex = ::floor(sTimer.elapsedTime() / interval);
2158
2159             if ((nIndex - oIndex) > VSMALL)
2160             {
2161                 oIndex = nIndex;
2162
2163                 scalar percent =
2164                 (
2165                     100.0 -
2166                     (
2167                         (100.0 * mesh.stack(tIndex).size())
2168                       / (stackSize + VSMALL)
2169                     )
2170                 );
2171
2172                 Info<< "  Swap Progress: " << percent << "% :"
2173                     << "  Total: " << mesh.status(1)
2174                     << "             \r"
2175                     << flush;
2176
2177                 reported = true;
2178             }
2179         }
2180
2181         // Retrieve the index for this face
2182         label fIndex = mesh.stack(tIndex).pop();
2183
2184         // Perform a Delaunay test and check if a flip is necesary.
2185         bool failed = mesh.testDelaunay(fIndex);
2186
2187         if (failed)
2188         {
2189             if (thread->master())
2190             {
2191                 // Swap this face.
2192                 mesh.swapQuadFace(fIndex);
2193             }
2194             else
2195             {
2196                 // Push this on to the master stack
2197                 mesh.stack(0).push(fIndex);
2198             }
2199         }
2200     }
2201
2202     if (thread->slave())
2203     {
2204         thread->sendSignal(meshHandler::STOP);
2205     }
2206
2207     if (reported)
2208     {
2209         Info<< "  Swap Progress: 100% :"
2210             << "  Total: " << mesh.status(1)
2211             << "             \r"
2212             << endl;
2213     }
2214 }
2215
2216
2217 // 3D Edge-swapping engine
2218 void dynamicTopoFvMesh::swap3DEdges
2219 (
2220     void *argument
2221 )
2222 {
2223     // Recast the argument
2224     meshHandler *thread = static_cast<meshHandler*>(argument);
2225
2226     if (thread->slave())
2227     {
2228         thread->sendSignal(meshHandler::START);
2229     }
2230
2231     dynamicTopoFvMesh& mesh = thread->reference();
2232
2233     // Figure out which thread this is...
2234     label tIndex = mesh.self();
2235
2236     // Dynamic programming variables
2237     labelList m;
2238     PtrList<scalarListList> Q;
2239     PtrList<labelListList> K, triangulations;
2240
2241     // Hull vertices information
2242     labelList hullV;
2243
2244     // Allocate dynamic programming tables
2245     mesh.initTables(m, Q, K, triangulations);
2246
2247     // Set the timer
2248     clockTime sTimer;
2249
2250     bool reported = false;
2251     label stackSize = mesh.stack(tIndex).size();
2252     scalar interval = mesh.reportInterval(), oIndex = 0.0, nIndex = 0.0;
2253
2254     oIndex = ::floor(sTimer.elapsedTime() / interval);
2255
2256     // Pick edges off the stack
2257     while (!mesh.stack(tIndex).empty())
2258     {
2259         // Report progress
2260         if (thread->master())
2261         {
2262             // Update the index, if its changed
2263             nIndex = ::floor(sTimer.elapsedTime() / interval);
2264
2265             if ((nIndex - oIndex) > VSMALL)
2266             {
2267                 oIndex = nIndex;
2268
2269                 scalar percent =
2270                 (
2271                     100.0 -
2272                     (
2273                         (100.0 * mesh.stack(tIndex).size())
2274                       / (stackSize + VSMALL)
2275                     )
2276                 );
2277
2278                 Info<< "  Swap Progress: " << percent << "% :"
2279                     << "  Surface: " << mesh.status(2)
2280                     << ", Total: " << mesh.status(1)
2281                     << "             \r"
2282                     << flush;
2283
2284                 reported = true;
2285             }
2286         }
2287
2288         // Retrieve an edge from the stack
2289         label eIndex = mesh.stack(tIndex).pop();
2290
2291         // Compute the minimum quality of cells around this edge
2292         scalar minQuality = mesh.computeMinQuality(eIndex, hullV);
2293
2294         // Check if this edge is on a bounding curve
2295         // (Override purity check for processor edges)
2296         if (mesh.checkBoundingCurve(eIndex, true))
2297         {
2298             continue;
2299         }
2300
2301         // Fill the dynamic programming tables
2302         if (mesh.fillTables(eIndex, minQuality, m, hullV, Q, K, triangulations))
2303         {
2304             // Check if edge-swapping is required.
2305             if (mesh.checkQuality(eIndex, m, Q, minQuality))
2306             {
2307                 if (thread->master())
2308                 {
2309                     // Remove this edge according to the swap sequence
2310                     mesh.removeEdgeFlips
2311                     (
2312                         eIndex,
2313                         minQuality,
2314                         hullV,
2315                         Q,
2316                         K,
2317                         triangulations
2318                     );
2319                 }
2320                 else
2321                 {
2322                     // Push this on to the master stack
2323                     mesh.stack(0).push(eIndex);
2324                 }
2325             }
2326         }
2327     }
2328
2329     if (thread->slave())
2330     {
2331         thread->sendSignal(meshHandler::STOP);
2332     }
2333
2334     if (reported)
2335     {
2336         Info<< "  Swap Progress: 100% :"
2337             << "  Surface: " << mesh.status(2)
2338             << ", Total: " << mesh.status(1)
2339             << "             \r"
2340             << endl;
2341     }
2342 }
2343
2344
2345 // Edge refinement engine
2346 void dynamicTopoFvMesh::edgeRefinementEngine
2347 (
2348     void *argument
2349 )
2350 {
2351     // Loop through all edges and bisect/collapse by the criterion:
2352     // Bisect when edge-length > ratioMax_*lengthScale
2353     // Collapse when edge-length < ratioMin_*lengthScale
2354
2355     // Recast the argument
2356     meshHandler *thread = static_cast<meshHandler*>(argument);
2357
2358     if (thread->slave())
2359     {
2360         thread->sendSignal(meshHandler::START);
2361     }
2362
2363     dynamicTopoFvMesh& mesh = thread->reference();
2364
2365     // Figure out which thread this is...
2366     label tIndex = mesh.self();
2367
2368     // Set the timer
2369     clockTime sTimer;
2370
2371     bool reported = false;
2372     label stackSize = mesh.stack(tIndex).size();
2373     scalar interval = mesh.reportInterval(), oIndex = 0.0, nIndex = 0.0;
2374
2375     oIndex = ::floor(sTimer.elapsedTime() / interval);
2376
2377     while (!mesh.stack(tIndex).empty())
2378     {
2379         // Update the index, if its changed
2380         // Report progress
2381         if (thread->master())
2382         {
2383             nIndex = ::floor(sTimer.elapsedTime() / interval);
2384
2385             if ((nIndex - oIndex) > VSMALL)
2386             {
2387                 oIndex = nIndex;
2388
2389                 scalar percent =
2390                 (
2391                     100.0 -
2392                     (
2393                         (100.0 * mesh.stack(tIndex).size())
2394                       / (stackSize + VSMALL)
2395                     )
2396                 );
2397
2398                 Info<< "  Refinement Progress: " << percent << "% :"
2399                     << "  Bisections: " << mesh.status(3)
2400                     << ", Collapses: " << mesh.status(4)
2401                     << ", Total: " << mesh.status(0)
2402                     << "             \r"
2403                     << flush;
2404
2405                 reported = true;
2406             }
2407         }
2408
2409         // Retrieve an entity from the stack
2410         label eIndex = mesh.stack(tIndex).pop();
2411
2412         if (mesh.checkBisection(eIndex))
2413         {
2414             if (thread->master())
2415             {
2416                 // Bisect this edge
2417                 mesh.bisectEdge(eIndex);
2418             }
2419             else
2420             {
2421                 // Push this on to the master stack
2422                 mesh.stack(0).push(eIndex);
2423             }
2424         }
2425         else
2426         if (mesh.checkCollapse(eIndex))
2427         {
2428             if (thread->master())
2429             {
2430                 // Collapse this edge
2431                 mesh.collapseEdge(eIndex);
2432             }
2433             else
2434             {
2435                 // Push this on to the master stack
2436                 mesh.stack(0).push(eIndex);
2437             }
2438         }
2439     }
2440
2441     if (thread->slave())
2442     {
2443         thread->sendSignal(meshHandler::STOP);
2444     }
2445
2446     if (reported)
2447     {
2448         Info<< "  Refinement Progress: 100% :"
2449             << "  Bisections: " << mesh.status(3)
2450             << ", Collapses: " << mesh.status(4)
2451             << ", Total: " << mesh.status(0)
2452             << "             \r"
2453             << endl;
2454     }
2455 }
2456
2457
2458 // Remove 2D sliver cells from the mesh
2459 void dynamicTopoFvMesh::remove2DSlivers()
2460 {
2461     // If coupled patches exist, set the flag
2462     if (patchCoupling_.size() || procIndices_.size())
2463     {
2464         setCoupledModification();
2465     }
2466
2467     // Sort by sliver-quality.
2468     labelList cIndices(thresholdSlivers_.toc());
2469     SortableList<scalar> values(cIndices.size());
2470
2471     // Fill-in values to sort by...
2472     forAll(cIndices, indexI)
2473     {
2474         values[indexI] = thresholdSlivers_[cIndices[indexI]];
2475     }
2476
2477     // Explicitly sort by quality value.
2478     values.sort();
2479
2480     const labelList& indices = values.indices();
2481
2482     if (debug)
2483     {
2484         if (thresholdSlivers_.size())
2485         {
2486             Pout<< "Sliver list: " << endl;
2487
2488             forAll(indices, indexI)
2489             {
2490                 label cIndex = cIndices[indices[indexI]];
2491
2492                 Pout<< " Cell: " << cIndex
2493                     << " Quality: " << thresholdSlivers_[cIndex]
2494                     << endl;
2495             }
2496
2497             if (debug > 1)
2498             {
2499                 writeVTK("sliverCells", cIndices, 3);
2500             }
2501         }
2502     }
2503
2504     forAll(indices, indexI)
2505     {
2506         // Fetch the cell
2507         label cIndex = cIndices[indices[indexI]];
2508
2509         // Ensure that this cell actually exists.
2510         if (cells_[cIndex].empty())
2511         {
2512             continue;
2513         }
2514
2515         const cell& cellToCheck = cells_[cIndex];
2516
2517         // Find an appropriate quad-face
2518         label fIndex = -1;
2519
2520         forAll(cellToCheck, faceI)
2521         {
2522             if (faces_[cellToCheck[faceI]].size() == 4)
2523             {
2524                 fIndex = cellToCheck[faceI];
2525                 break;
2526             }
2527         }
2528
2529         // Find the prism faces
2530         FixedList<label,2> cBdyIndex(-1), cIntIndex(-1);
2531         FixedList<face,2> cBdyFace, cIntFace;
2532
2533         meshOps::findPrismFaces
2534         (
2535             fIndex,
2536             cIndex,
2537             faces_,
2538             cells_,
2539             neighbour_,
2540             cBdyFace,
2541             cBdyIndex,
2542             cIntFace,
2543             cIntIndex
2544         );
2545
2546         if (debug > 1)
2547         {
2548             InfoIn("void dynamicTopoFvMesh::remove2DSlivers()")
2549                 << nl
2550                 << " Considering cell: " << cIndex
2551                 << ":: " << cells_[cIndex]
2552                 << " for sliver removal."
2553                 << " Using face: " << fIndex
2554                 << ":: " << faces_[fIndex]
2555                 << endl;
2556         }
2557
2558         label triFace = cBdyIndex[0], triEdge = -1;
2559
2560         // Find the common-edge between quad-tri faces
2561         meshOps::findCommonEdge
2562         (
2563             fIndex,
2564             triFace,
2565             faceEdges_,
2566             triEdge
2567         );
2568
2569         // Find the isolated point.
2570         label ptIndex = -1, nextPtIndex = -1;
2571
2572         const edge& edgeToCheck = edges_[triEdge];
2573
2574         meshOps::findIsolatedPoint
2575         (
2576             faces_[triFace],
2577             edgeToCheck,
2578             ptIndex,
2579             nextPtIndex
2580         );
2581
2582         // Determine the interior faces connected to each edge-point.
2583         label firstFace = -1, secondFace = -1;
2584
2585         if (cIntFace[0].which(edgeToCheck[0]) > -1)
2586         {
2587             firstFace  = cIntIndex[0];
2588             secondFace = cIntIndex[1];
2589         }
2590         else
2591         {
2592             firstFace  = cIntIndex[1];
2593             secondFace = cIntIndex[0];
2594         }
2595
2596         point ec =
2597         (
2598             linePointRef
2599             (
2600                 points_[edgeToCheck[0]],
2601                 points_[edgeToCheck[1]]
2602             ).centre()
2603         );
2604
2605         FixedList<vector, 2> p(vector::zero), q(vector::zero);
2606         FixedList<scalar, 2> proj(0.0);
2607
2608         // Find the projection on the edge.
2609         forAll(edgeToCheck, pointI)
2610         {
2611             p[pointI] = (points_[ptIndex] - points_[edgeToCheck[pointI]]);
2612             q[pointI] = (ec - points_[edgeToCheck[pointI]]);
2613
2614             q[pointI] /= (mag(q[pointI]) + VSMALL);
2615
2616             proj[pointI] = (p[pointI] & q[pointI]);
2617         }
2618
2619         // Take action based on the magnitude of the projection.
2620         if (mag(proj[0]) < VSMALL)
2621         {
2622             if (collapseQuadFace(firstFace).type() > 0)
2623             {
2624                 statistics_[7]++;
2625             }
2626
2627             continue;
2628         }
2629
2630         if (mag(proj[1]) < VSMALL)
2631         {
2632             if (collapseQuadFace(secondFace).type() > 0)
2633             {
2634                 statistics_[7]++;
2635             }
2636
2637             continue;
2638         }
2639
2640         if (proj[0] > 0.0 && proj[1] < 0.0)
2641         {
2642             changeMap map = bisectQuadFace(firstFace, changeMap());
2643
2644             // Loop through added faces, and collapse
2645             // the appropriate one
2646             const List<objectMap>& aF = map.addedFaceList();
2647
2648             forAll(aF, faceI)
2649             {
2650                 if
2651                 (
2652                     (owner_[aF[faceI].index()] == cIndex) &&
2653                     (aF[faceI].index() != firstFace)
2654                 )
2655                 {
2656                     if (collapseQuadFace(aF[faceI].index()).type() > 0)
2657                     {
2658                         statistics_[7]++;
2659                     }
2660
2661                     break;
2662                 }
2663             }
2664
2665             continue;
2666         }
2667
2668         if (proj[0] < 0.0 && proj[1] > 0.0)
2669         {
2670             changeMap map = bisectQuadFace(secondFace, changeMap());
2671
2672             // Loop through added faces, and collapse
2673             // the appropriate one
2674             const List<objectMap>& aF = map.addedFaceList();
2675
2676             forAll(aF, faceI)
2677             {
2678                 if
2679                 (
2680                     (owner_[aF[faceI].index()] == cIndex) &&
2681                     (aF[faceI].index() != secondFace)
2682                 )
2683                 {
2684                     if (collapseQuadFace(aF[faceI].index()).type() > 0)
2685                     {
2686                         statistics_[7]++;
2687                     }
2688
2689                     break;
2690                 }
2691             }
2692
2693             continue;
2694         }
2695
2696         if (proj[0] > 0.0 && proj[1] > 0.0)
2697         {
2698             changeMap map = bisectQuadFace(fIndex, changeMap());
2699
2700             // Loop through added faces, and collapse
2701             // the appropriate one
2702             const List<objectMap>& aF = map.addedFaceList();
2703
2704             forAll(aF, faceI)
2705             {
2706                 if
2707                 (
2708                     (owner_[aF[faceI].index()] == cIndex) &&
2709                     (aF[faceI].index() != fIndex)
2710                 )
2711                 {
2712                     if (collapseQuadFace(aF[faceI].index()).type() > 0)
2713                     {
2714                         statistics_[7]++;
2715                     }
2716
2717                     break;
2718                 }
2719             }
2720
2721             continue;
2722         }
2723     }
2724
2725     // Clear out the list
2726     thresholdSlivers_.clear();
2727
2728     // If coupled patches exist, reset the flag
2729     if (patchCoupling_.size() || procIndices_.size())
2730     {
2731         unsetCoupledModification();
2732     }
2733 }
2734
2735
2736 // Identify the sliver type in 3D
2737 const changeMap dynamicTopoFvMesh::identifySliverType
2738 (
2739     const label cIndex
2740 ) const
2741 {
2742     changeMap map;
2743
2744     // Ensure that this cell actually exists.
2745     if (cells_[cIndex].empty())
2746     {
2747         return map;
2748     }
2749
2750     label fourthPoint = -1;
2751     scalar minDistance = GREAT;
2752     face tFace(3), testFace(3), faceToCheck(3);
2753     FixedList<edge, 2> edgeToCheck(edge(-1,-1));
2754
2755     const cell& cellToCheck = cells_[cIndex];
2756
2757     // Find the point-face pair with minimum perpendicular distance
2758     forAll(cellToCheck, faceI)
2759     {
2760         label isolatedPoint = -1;
2761         label nextFaceI = cellToCheck.fcIndex(faceI);
2762
2763         // Pick two faces from this cell.
2764         const face& currFace = faces_[cellToCheck[faceI]];
2765         const face& nextFace = faces_[cellToCheck[nextFaceI]];
2766
2767         // Get the fourth point
2768         forAll(nextFace, pointI)
2769         {
2770             if
2771             (
2772                 nextFace[pointI] != currFace[0]
2773              && nextFace[pointI] != currFace[1]
2774              && nextFace[pointI] != currFace[2]
2775             )
2776             {
2777                 isolatedPoint = nextFace[pointI];
2778
2779                 // Configure a triangular face with correct orientation.
2780                 if (owner_[cellToCheck[faceI]] == cIndex)
2781                 {
2782                     testFace[0] = currFace[2];
2783                     testFace[1] = currFace[1];
2784                     testFace[2] = currFace[0];
2785                 }
2786                 else
2787                 {
2788                     testFace[0] = currFace[0];
2789                     testFace[1] = currFace[1];
2790                     testFace[2] = currFace[2];
2791                 }
2792
2793                 break;
2794             }
2795         }
2796
2797         // Obtain the unit normal.
2798         vector testNormal = testFace.normal(points_);
2799
2800         testNormal /= (mag(testNormal) + VSMALL);
2801
2802         // Project the isolated point onto the face.
2803         vector p = points_[isolatedPoint] - points_[testFace[0]];
2804         vector q = p - ((p & testNormal)*testNormal);
2805
2806         // Compute the distance
2807         scalar distance = mag(p - q);
2808
2809         // Is it the least so far?
2810         if (distance < minDistance)
2811         {
2812             // Use this point-face pair.
2813             fourthPoint = isolatedPoint;
2814             tFace = testFace;
2815             minDistance = distance;
2816         }
2817     }
2818
2819     // Obtain the face-normal.
2820     vector refArea = tFace.normal(points_);
2821
2822     // Normalize it.
2823     vector n = refArea/mag(refArea);
2824
2825     // Define edge-vectors.
2826     vector r1 = points_[tFace[1]] - points_[tFace[0]];
2827     vector r2 = points_[tFace[2]] - points_[tFace[1]];
2828     vector r3 = points_[tFace[0]] - points_[tFace[2]];
2829
2830     // Project the fourth point onto the face.
2831     vector r4 = points_[fourthPoint] - points_[tFace[0]];
2832
2833     r4 = r4 - ((r4 & n)*n);
2834
2835     // Define the two other vectors.
2836     vector r5 = r4 - r1;
2837     vector r6 = r5 - r2;
2838
2839     // Calculate three signed triangle areas, using tFace[0] as the origin.
2840     scalar t1 = n & (0.5 * (r1 ^ r4));
2841     scalar t2 = n & (0.5 * (r2 ^ r5));
2842     scalar t3 = n & (0.5 * (r3 ^ r6));
2843
2844     // Determine sliver types based on are magnitudes.
2845     if (t1 > 0 && t2 > 0 && t3 > 0)
2846     {
2847         // Region R0: Cap cell.
2848         map.type() = 2;
2849         map.add("apexPoint", fourthPoint);
2850
2851         faceToCheck[0] = tFace[0];
2852         faceToCheck[1] = tFace[1];
2853         faceToCheck[2] = tFace[2];
2854     }
2855
2856     if (t1 < 0 && t2 > 0 && t3 > 0)
2857     {
2858         // Region R1: Sliver cell.
2859         map.type() = 1;
2860
2861         edgeToCheck[0][0] = tFace[2];
2862         edgeToCheck[0][1] = fourthPoint;
2863         edgeToCheck[1][0] = tFace[0];
2864         edgeToCheck[1][1] = tFace[1];
2865     }
2866
2867     if (t1 > 0 && t2 < 0 && t3 > 0)
2868     {
2869         // Region R2: Sliver cell.
2870         map.type() = 1;
2871
2872         edgeToCheck[0][0] = tFace[0];
2873         edgeToCheck[0][1] = fourthPoint;
2874         edgeToCheck[1][0] = tFace[1];
2875         edgeToCheck[1][1] = tFace[2];
2876     }
2877
2878     if (t1 > 0 && t2 > 0 && t3 < 0)
2879     {
2880         // Region R3: Sliver cell.
2881         map.type() = 1;
2882
2883         edgeToCheck[0][0] = tFace[1];
2884         edgeToCheck[0][1] = fourthPoint;
2885         edgeToCheck[1][0] = tFace[2];
2886         edgeToCheck[1][1] = tFace[0];
2887     }
2888
2889     if (t1 < 0 && t2 > 0 && t3 < 0)
2890     {
2891         // Region R4: Cap cell.
2892         map.type() = 2;
2893         map.add("apexPoint", tFace[0]);
2894
2895         faceToCheck[0] = tFace[1];
2896         faceToCheck[1] = tFace[2];
2897         faceToCheck[2] = fourthPoint;
2898     }
2899
2900     if (t1 < 0 && t2 < 0 && t3 > 0)
2901     {
2902         // Region R5: Cap cell.
2903         map.type() = 2;
2904         map.add("apexPoint", tFace[1]);
2905
2906         faceToCheck[0] = tFace[2];
2907         faceToCheck[1] = tFace[0];
2908         faceToCheck[2] = fourthPoint;
2909     }
2910
2911     if (t1 > 0 && t2 < 0 && t3 < 0)
2912     {
2913         // Region R6: Cap cell.
2914         map.type() = 2;
2915         map.add("apexPoint", tFace[2]);
2916
2917         faceToCheck[0] = tFace[0];
2918         faceToCheck[1] = tFace[1];
2919         faceToCheck[2] = fourthPoint;
2920     }
2921
2922     // See if an over-ride to wedge/spade is necessary.
2923     // Obtain a reference area magnitude.
2924     scalar refMag = 0.1*(refArea & n);
2925
2926     if (mag(t1) < refMag)
2927     {
2928         if (mag(t3) < refMag)
2929         {
2930             // Wedge case: Too close to point [0]
2931             map.type() = 4;
2932
2933             edgeToCheck[0][0] = tFace[0];
2934             edgeToCheck[0][1] = fourthPoint;
2935         }
2936         else
2937         if (mag(t2) < refMag)
2938         {
2939             // Wedge case: Too close to point [1]
2940             map.type() = 4;
2941
2942             edgeToCheck[0][0] = tFace[1];
2943             edgeToCheck[0][1] = fourthPoint;
2944         }
2945         else
2946         if ((mag(t2) > refMag) && (mag(t3) > refMag))
2947         {
2948             // Spade case: Too close to edge vector r1
2949             map.type() = 3;
2950             map.add("apexPoint", fourthPoint);
2951
2952             edgeToCheck[0][0] = tFace[0];
2953             edgeToCheck[0][1] = tFace[1];
2954         }
2955     }
2956
2957     if (mag(t2) < refMag)
2958     {
2959         if (mag(t3) < refMag)
2960         {
2961             // Wedge case: Too close to point [2]
2962             map.type() = 4;
2963
2964             edgeToCheck[0][0] = tFace[2];
2965             edgeToCheck[0][1] = fourthPoint;
2966         }
2967         else
2968         if ((mag(t1) > refMag) && (mag(t3) > refMag))
2969         {
2970             // Spade case: Too close to edge vector r2
2971             map.type() = 3;
2972             map.add("apexPoint", fourthPoint);
2973
2974             edgeToCheck[0][0] = tFace[1];
2975             edgeToCheck[0][1] = tFace[2];
2976         }
2977     }
2978
2979     if (mag(t3) < refMag)
2980     {
2981         if ((mag(t1) > refMag) && (mag(t2) > refMag))
2982         {
2983             // Spade case: Too close to edge vector r3
2984             map.type() = 3;
2985             map.add("apexPoint", fourthPoint);
2986
2987             edgeToCheck[0][0] = tFace[2];
2988             edgeToCheck[0][1] = tFace[0];
2989         }
2990     }
2991
2992     // Determine appropriate information for sliver exudation.
2993     switch (map.type())
2994     {
2995         case 1:
2996         {
2997             FixedList<bool, 2> foundEdge(false);
2998
2999             // Search the cell-faces for first and second edges.
3000             forAll(cellToCheck, faceI)
3001             {
3002                 const labelList& fEdges = faceEdges_[cellToCheck[faceI]];
3003
3004                 forAll(fEdges, edgeI)
3005                 {
3006                     const edge& thisEdge = edges_[fEdges[edgeI]];
3007
3008                     if (thisEdge == edgeToCheck[0])
3009                     {
3010                         map.add("firstEdge", fEdges[edgeI]);
3011
3012                         foundEdge[0] = true;
3013                     }
3014
3015                     if (thisEdge == edgeToCheck[1])
3016                     {
3017                         map.add("secondEdge", fEdges[edgeI]);
3018
3019                         foundEdge[1] = true;
3020                     }
3021                 }
3022
3023                 if (foundEdge[0] && foundEdge[1])
3024                 {
3025                     break;
3026                 }
3027             }
3028
3029             break;
3030         }
3031
3032         case 2:
3033         {
3034             // Search the cell-faces for opposing faces.
3035             forAll(cellToCheck, faceI)
3036             {
3037                 const face& thisFace = faces_[cellToCheck[faceI]];
3038
3039                 if (triFace::compare(triFace(thisFace), triFace(faceToCheck)))
3040                 {
3041                     map.add("opposingFace", cellToCheck[faceI]);
3042
3043                     break;
3044                 }
3045             }
3046
3047             break;
3048         }
3049
3050         case 3:
3051         case 4:
3052         {
3053             bool foundEdge = false;
3054
3055             // Search the cell-faces for first edge.
3056             forAll(cellToCheck, faceI)
3057             {
3058                 const labelList& fEdges = faceEdges_[cellToCheck[faceI]];
3059
3060                 forAll(fEdges, edgeI)
3061                 {
3062                     const edge& thisEdge = edges_[fEdges[edgeI]];
3063
3064                     if (thisEdge == edgeToCheck[0])
3065                     {
3066                         map.add("firstEdge", fEdges[edgeI]);
3067
3068                         foundEdge = true;
3069                     }
3070                 }
3071
3072                 if (foundEdge)
3073                 {
3074                     break;
3075                 }
3076             }
3077
3078             break;
3079         }
3080
3081         default:
3082         {
3083             WarningIn
3084             (
3085                 "void dynamicTopoFvMesh::identifySliverType"
3086                 "(const label cIndex) const"
3087             )
3088                 << nl << "Could not identify sliver type for cell: "
3089                 << cIndex
3090                 << endl;
3091         }
3092     }
3093
3094     if (debug > 2)
3095     {
3096         Pout<< "Cell: " << cIndex
3097             << " Identified sliver type as: "
3098             << map.type() << endl;
3099     }
3100
3101     // Return the result.
3102     return map;
3103 }
3104
3105
3106 // Remove sliver cells
3107 void dynamicTopoFvMesh::removeSlivers()
3108 {
3109     if (!edgeRefinement_)
3110     {
3111         return;
3112     }
3113
3114     // Check if a removeSlivers entry was found in the dictionary
3115     if (dict_.subDict("dynamicTopoFvMesh").found("removeSlivers"))
3116     {
3117         Switch rs =
3118         (
3119             dict_.subDict("dynamicTopoFvMesh").lookup("removeSlivers")
3120         );
3121
3122         if (!rs)
3123         {
3124             return;
3125         }
3126     }
3127
3128     // Invoke the 2D sliver removal routine
3129     if (twoDMesh_)
3130     {
3131         remove2DSlivers();
3132         return;
3133     }
3134
3135     // If coupled patches exist, set the flag
3136     if (patchCoupling_.size() || procIndices_.size())
3137     {
3138         setCoupledModification();
3139     }
3140
3141     // Sort by sliver-quality.
3142     labelList cIndices(thresholdSlivers_.toc());
3143     SortableList<scalar> values(cIndices.size());
3144
3145     // Fill-in values to sort by...
3146     forAll(cIndices, indexI)
3147     {
3148         values[indexI] = thresholdSlivers_[cIndices[indexI]];
3149     }
3150
3151     // Explicitly sort by quality value.
3152     values.sort();
3153
3154     const labelList& indices = values.indices();
3155
3156     if (debug && thresholdSlivers_.size())
3157     {
3158         Pout<< "Sliver list: " << endl;
3159
3160         forAll(indices, indexI)
3161         {
3162             label cIndex = cIndices[indices[indexI]];
3163
3164             Pout<< " Cell: " << cIndex
3165                 << " Quality: " << thresholdSlivers_[cIndex]
3166                 << endl;
3167         }
3168
3169         if (debug > 1)
3170         {
3171             writeVTK("sliverCells", cIndices, 3);
3172         }
3173     }
3174
3175     forAll(indices, indexI)
3176     {
3177         // Fetch the cell index
3178         label cIndex = cIndices[indices[indexI]];
3179
3180         // First check if this sliver cell is handled elsewhere.
3181         if (procIndices_.size())
3182         {
3183             bool foundInSubMesh = false;
3184
3185             forAll(procIndices_, procI)
3186             {
3187                 label proc = procIndices_[procI];
3188
3189                 if (proc < Pstream::myProcNo())
3190                 {
3191                     Map<label>& rCellMap =
3192                     (
3193                         sendMeshes_[proc].map().reverseEntityMap
3194                         (
3195                             coupleMap::CELL
3196                         )
3197                     );
3198
3199                     if (rCellMap.found(cIndex))
3200                     {
3201                         // This cell was sent to another sub-domain.
3202                         foundInSubMesh = true;
3203                         break;
3204                     }
3205                 }
3206             }
3207
3208             if (foundInSubMesh)
3209             {
3210                 continue;
3211             }
3212         }
3213
3214         // Identify the sliver type.
3215         changeMap map = identifySliverType(cIndex);
3216
3217         if (debug)
3218         {
3219             InfoIn("void dynamicTopoFvMesh::removeSlivers()")
3220                 << nl << "Removing Cell: " << cIndex
3221                 << " of sliver type: " << map.type()
3222                 << " with quality: " << thresholdSlivers_[cIndex]
3223                 << endl;
3224         }
3225
3226         bool success = false;
3227
3228         // Take action based on the type of sliver.
3229         switch (map.type())
3230         {
3231             case -1:
3232             {
3233                 // Sliver cell was removed by a prior operation.
3234                 // Nothing needs to be done.
3235                 break;
3236             }
3237
3238             case 1:
3239             {
3240                 // Sliver cell.
3241                 // Determine which edges need to be bisected.
3242                 label firstEdge = readLabel(map.lookup("firstEdge"));
3243                 label secondEdge = readLabel(map.lookup("secondEdge"));
3244
3245                 // Force bisection on both edges.
3246                 changeMap firstMap  = bisectEdge(firstEdge, false, true);
3247                 changeMap secondMap = bisectEdge(secondEdge, false, true);
3248
3249                 // Collapse the intermediate edge.
3250                 // Since we don't know which edge it is, search
3251                 // through recently added edges and compare.
3252                 edge edgeToCheck
3253                 (
3254                     firstMap.addedPointList()[0].index(),
3255                     secondMap.addedPointList()[0].index()
3256                 );
3257
3258                 bool foundCollapseEdge = false;
3259
3260                 const List<objectMap>& firstMapEdges =
3261                 (
3262                     firstMap.addedEdgeList()
3263                 );
3264
3265                 const List<objectMap>& secondMapEdges =
3266                 (
3267                     secondMap.addedEdgeList()
3268                 );
3269
3270                 // Loop through the first list.
3271                 forAll(firstMapEdges, edgeI)
3272                 {
3273                     const edge& thisEdge =
3274                     (
3275                         edges_[firstMapEdges[edgeI].index()]
3276                     );
3277
3278                     if (thisEdge == edgeToCheck)
3279                     {
3280                         // Collapse this edge.
3281                         if
3282                         (
3283                             collapseEdge
3284                             (
3285                                 firstMapEdges[edgeI].index(),
3286                                 -1,
3287                                 false,
3288                                 true
3289                             ).type() > 0
3290                         )
3291                         {
3292                             success = true;
3293                         }
3294
3295                         foundCollapseEdge = true;
3296                         break;
3297                     }
3298                 }
3299
3300                 // Loop through the second list.
3301                 if (!foundCollapseEdge)
3302                 {
3303                     forAll(secondMapEdges, edgeI)
3304                     {
3305                         const edge& thisEdge =
3306                         (
3307                             edges_[secondMapEdges[edgeI].index()]
3308                         );
3309
3310                         if (thisEdge == edgeToCheck)
3311                         {
3312                             // Collapse this edge.
3313                             if
3314                             (
3315                                 collapseEdge
3316                                 (
3317                                     secondMapEdges[edgeI].index(),
3318                                     -1,
3319                                     false,
3320                                     true
3321                                 ).type() > 0
3322                             )
3323                             {
3324                                 success = true;
3325                             }
3326
3327                             break;
3328                         }
3329                     }
3330                 }
3331
3332                 break;
3333             }
3334
3335             case 2:
3336             {
3337                 // Cap cell.
3338                 label opposingFace = readLabel(map.lookup("opposingFace"));
3339
3340                 // Force trisection of the opposing face.
3341                 changeMap faceMap =
3342                 (
3343                     trisectFace(opposingFace, false, true)
3344                 );
3345
3346                 // Collapse the intermediate edge.
3347                 // Since we don't know which edge it is, search
3348                 // through recently added edges and compare.
3349                 edge edgeToCheck
3350                 (
3351                     readLabel(map.lookup("apexPoint")),
3352                     faceMap.addedPointList()[0].index()
3353                 );
3354
3355                 const List<objectMap>& faceMapEdges =
3356                 (
3357                     faceMap.addedEdgeList()
3358                 );
3359
3360                 forAll(faceMapEdges, edgeI)
3361                 {
3362                     const edge& thisEdge = edges_[faceMapEdges[edgeI].index()];
3363
3364                     if (thisEdge == edgeToCheck)
3365                     {
3366                         // Collapse this edge.
3367                         if
3368                         (
3369                             collapseEdge
3370                             (
3371                                 faceMapEdges[edgeI].index(),
3372                                 -1,
3373                                 false,
3374                                 true
3375                             ).type() > 0
3376                         )
3377                         {
3378                             success = true;
3379                         }
3380
3381                         break;
3382                     }
3383                 }
3384
3385                 break;
3386             }
3387
3388             case 3:
3389             {
3390                 // Spade cell.
3391
3392                 // Force bisection on the first edge.
3393                 changeMap firstMap =
3394                 (
3395                     bisectEdge
3396                     (
3397                         readLabel(map.lookup("firstEdge")),
3398                         false,
3399                         true
3400                     )
3401                 );
3402
3403                 // Collapse the intermediate edge.
3404                 // Since we don't know which edge it is, search
3405                 // through recently added edges and compare.
3406                 edge edgeToCheck
3407                 (
3408                     readLabel(map.lookup("apexPoint")),
3409                     firstMap.addedPointList()[0].index()
3410                 );
3411
3412                 const List<objectMap>& firstMapEdges =
3413                 (
3414                     firstMap.addedEdgeList()
3415                 );
3416
3417                 // Loop through the first list.
3418                 forAll(firstMapEdges, edgeI)
3419                 {
3420                     const edge& thisEdge = edges_[firstMapEdges[edgeI].index()];
3421
3422                     if (thisEdge == edgeToCheck)
3423                     {
3424                         // Collapse this edge.
3425                         if
3426                         (
3427                             collapseEdge
3428                             (
3429                                 firstMapEdges[edgeI].index(),
3430                                 -1,
3431                                 false,
3432                                 true
3433                             ).type() > 0
3434                         )
3435                         {
3436                             success = true;
3437                         }
3438
3439                         break;
3440                     }
3441                 }
3442
3443                 break;
3444             }
3445
3446             case 4:
3447             {
3448                 // Wedge cell.
3449
3450                 // Collapse the first edge.
3451                 if
3452                 (
3453                     collapseEdge
3454                     (
3455                         readLabel(map.lookup("firstEdge")),
3456                         -1,
3457                         false,
3458                         true
3459                     ).type() > 0
3460                 )
3461                 {
3462                     success = true;
3463                 }
3464
3465                 break;
3466             }
3467
3468             default:
3469             {
3470                 WarningIn("void dynamicTopoFvMesh::removeSlivers()")
3471                     << nl << "Could not identify sliver type for cell: "
3472                     << cIndex
3473                     << endl;
3474             }
3475         }
3476
3477         // Increment the count for successful sliver removal
3478         if (success)
3479         {
3480             statistics_[7]++;
3481         }
3482     }
3483
3484     // Clear out the list
3485     thresholdSlivers_.clear();
3486
3487     // If coupled patches exist, reset the flag
3488     if (patchCoupling_.size() || procIndices_.size())
3489     {
3490         unsetCoupledModification();
3491     }
3492 }
3493
3494
3495 // Given an input quad-face, determine checkEdges from mesh
3496 //  - Does not refer to member data directly,
3497 //    since this is also used by subMeshes
3498 void dynamicTopoFvMesh::getCheckEdges
3499 (
3500     const label fIndex,
3501     const dynamicTopoFvMesh& mesh,
3502     changeMap& map,
3503     FixedList<edge,4>& checkEdge,
3504     FixedList<label,4>& checkEdgeIndex
3505 )
3506 {
3507     checkEdgeIndex[0] = mesh.getTriBoundaryEdge(fIndex);
3508     checkEdge[0] = mesh.edges_[checkEdgeIndex[0]];
3509
3510     const labelList& fEdges = mesh.faceEdges_[fIndex];
3511
3512     forAll(fEdges, edgeI)
3513     {
3514         if (checkEdgeIndex[0] != fEdges[edgeI])
3515         {
3516             const edge& thisEdge = mesh.edges_[fEdges[edgeI]];
3517
3518             if
3519             (
3520                 checkEdge[0].start() == thisEdge[0] ||
3521                 checkEdge[0].start() == thisEdge[1]
3522             )
3523             {
3524                 checkEdgeIndex[1] = fEdges[edgeI];
3525                 checkEdge[1] = thisEdge;
3526
3527                 // Update the map
3528                 map.add("firstEdge", checkEdgeIndex[1]);
3529             }
3530             else
3531             if
3532             (
3533                 checkEdge[0].end() == thisEdge[0] ||
3534                 checkEdge[0].end() == thisEdge[1]
3535             )
3536             {
3537                 checkEdgeIndex[2] = fEdges[edgeI];
3538                 checkEdge[2] = thisEdge;
3539
3540                 // Update the map
3541                 map.add("secondEdge", checkEdgeIndex[2]);
3542             }
3543             else
3544             {
3545                 checkEdgeIndex[3] = fEdges[edgeI];
3546                 checkEdge[3] = thisEdge;
3547             }
3548         }
3549     }
3550 }
3551
3552
3553 // Return length-scale at an face-location in the mesh [2D]
3554 scalar dynamicTopoFvMesh::faceLengthScale
3555 (
3556     const label fIndex
3557 ) const
3558 {
3559     // Reset the scale first
3560     scalar scale = 0.0;
3561
3562     label facePatch = whichPatch(fIndex);
3563
3564     // Determine whether the face is internal
3565     if (facePatch < 0)
3566     {
3567         scale =
3568         (
3569             0.5 *
3570             (
3571                 lengthScale_[owner_[fIndex]]
3572               + lengthScale_[neighbour_[fIndex]]
3573             )
3574         );
3575     }
3576     else
3577     {
3578         // Fetch the fixed-length scale
3579         scale = lengthEstimator().fixedLengthScale(fIndex, facePatch);
3580
3581         // If this is a floating face, pick the owner length-scale
3582         if (lengthEstimator().isFreePatch(facePatch))
3583         {
3584             scale = lengthScale_[owner_[fIndex]];
3585         }
3586
3587         // If proximity-based refinement is requested,
3588         // test the proximity to the nearest non-neighbour.
3589         if (lengthEstimator().isProximityPatch(facePatch))
3590         {
3591             label proximityFace = -1;
3592
3593             // Perform a proximity-check.
3594             scalar distance = testProximity(fIndex, proximityFace);
3595
3596             if (debug > 3 && self() == 0)
3597             {
3598                 if
3599                 (
3600                     (proximityFace > -1) &&
3601                     ((distance / 5.0) < scale)
3602                 )
3603                 {
3604                     Pout<< " Closest opposing face detected for face: " << nl
3605                         << '\t' << fIndex
3606                         << " :: " << faces_[fIndex]
3607                         << " was face:\n"
3608                         << '\t' << proximityFace
3609                         << " :: " << polyMesh::faces()[proximityFace] << nl
3610                         << " with distance: " << distance
3611                         << endl;
3612                 }
3613             }
3614
3615             scale =
3616             (
3617                 Foam::min
3618                 (
3619                     scale,
3620                     ((distance / 3.0) - SMALL)/lengthEstimator().ratioMax()
3621                 )
3622             );
3623         }
3624
3625         // If this face lies on a processor patch,
3626         // fetch lengthScale info from patchSubMeshes
3627         if (processorCoupledEntity(fIndex))
3628         {
3629             scale = processorLengthScale(fIndex);
3630         }
3631
3632         // Limit scales if necessary
3633         lengthEstimator().limitScale(scale);
3634     }
3635
3636     return scale;
3637 }
3638
3639
3640 // Compute length-scale at an edge-location in the mesh [3D]
3641 scalar dynamicTopoFvMesh::edgeLengthScale
3642 (
3643     const label eIndex
3644 ) const
3645 {
3646     // Reset the scale first
3647     scalar scale = 0.0;
3648
3649     const labelList& eFaces = edgeFaces_[eIndex];
3650
3651     label edgePatch = whichEdgePatch(eIndex);
3652
3653     // Determine whether the edge is internal
3654     if (edgePatch < 0)
3655     {
3656         forAll(eFaces, faceI)
3657         {
3658             scale += lengthScale_[owner_[eFaces[faceI]]];
3659             scale += lengthScale_[neighbour_[eFaces[faceI]]];
3660         }
3661
3662         scale /= (2.0*eFaces.size());
3663     }
3664     else
3665     {
3666         // Search for boundary faces, and average their scale
3667         forAll(eFaces, faceI)
3668         {
3669             label facePatch = whichPatch(eFaces[faceI]);
3670
3671             // Skip internal faces
3672             if (facePatch == -1)
3673             {
3674                 continue;
3675             }
3676
3677             // If this is a floating face, pick the owner length-scale
3678             if (lengthEstimator().isFreePatch(facePatch))
3679             {
3680                 scale += lengthScale_[owner_[eFaces[faceI]]];
3681             }
3682             else
3683             {
3684                 // Fetch fixed length-scale
3685                 scale +=
3686                 (
3687                     lengthEstimator().fixedLengthScale
3688                     (
3689                         eFaces[faceI],
3690                         facePatch
3691                     )
3692                 );
3693             }
3694         }
3695
3696         scale *= 0.5;
3697
3698         // If proximity-based refinement is requested,
3699         // test the proximity to the nearest non-neighbour.
3700         if (lengthEstimator().isProximityPatch(edgePatch))
3701         {
3702             label proximityFace = -1;
3703
3704             // Perform a proximity-check.
3705             scalar distance = testProximity(eIndex, proximityFace);
3706
3707             if (debug > 3 && self() == 0)
3708             {
3709                 if
3710                 (
3711                     (proximityFace > -1) &&
3712                     ((distance / 5.0) < scale)
3713                 )
3714                 {
3715                     Pout<< " Closest opposing face detected for edge: " << nl
3716                         << '\t' << eIndex
3717                         << " :: " << edges_[eIndex]
3718                         << " was face:\n"
3719                         << '\t' << proximityFace
3720                         << " :: " << polyMesh::faces()[proximityFace] << nl
3721                         << " with distance: " << distance
3722                         << endl;
3723                 }
3724             }
3725
3726             scale =
3727             (
3728                 Foam::min
3729                 (
3730                     scale,
3731                     ((distance / 3.0) - SMALL)/lengthEstimator().ratioMax()
3732                 )
3733             );
3734         }
3735
3736         // If curvature-based refinement is requested,
3737         // test the variation in face-normal directions.
3738         if (lengthEstimator().isCurvaturePatch(edgePatch))
3739         {
3740             // Obtain face-normals for both faces.
3741             label count = 0;
3742             FixedList<vector, 2> fNorm;
3743
3744             forAll(eFaces, faceI)
3745             {
3746                 if (neighbour_[eFaces[faceI]] == -1)
3747                 {
3748                     // Obtain the normal.
3749                     fNorm[count] = faces_[eFaces[faceI]].normal(points_);
3750
3751                     // Normalize it.
3752                     fNorm[count] /= mag(fNorm[count]);
3753
3754                     count++;
3755                 }
3756             }
3757
3758             scalar deviation = (fNorm[0] & fNorm[1]);
3759             scalar refDeviation = lengthEstimator().curvatureDeviation();
3760
3761             if (mag(deviation) < refDeviation)
3762             {
3763                 // Fetch the edge
3764                 const edge& edgeToCheck = edges_[eIndex];
3765
3766                 // Get the edge-length.
3767                 scalar length =
3768                 (
3769                     linePointRef
3770                     (
3771                         points_[edgeToCheck.start()],
3772                         points_[edgeToCheck.end()]
3773                     ).mag()
3774                 );
3775
3776                 if (debug > 3 && self() == 0)
3777                 {
3778                     Pout<< "Deviation: " << deviation << nl
3779                         << "curvatureDeviation: " << refDeviation
3780                         << ", Edge: " << eIndex << ", Length: " << length
3781                         << ", Scale: " << scale << nl
3782                         << " Half-length: " << (0.5*length) << nl
3783                         << " MinRatio: "
3784                         << (lengthEstimator().ratioMin()*scale)
3785                         << endl;
3786                 }
3787
3788                 scale =
3789                 (
3790                     Foam::min
3791                     (
3792                         scale,
3793                         ((length - SMALL)/lengthEstimator().ratioMax())
3794                     )
3795                 );
3796             }
3797         }
3798
3799         // If this edge lies on a processor patch,
3800         // fetch lengthScale info from patchSubMeshes
3801         if (processorCoupledEntity(eIndex))
3802         {
3803             scale = processorLengthScale(eIndex);
3804         }
3805
3806         // Limit scales if necessary
3807         lengthEstimator().limitScale(scale);
3808     }
3809
3810     return scale;
3811 }
3812
3813
3814 // MultiThreaded topology modifier
3815 void dynamicTopoFvMesh::threadedTopoModifier()
3816 {
3817     // Remove sliver cells first.
3818     removeSlivers();
3819
3820     // Coupled entities to avoid during normal modification
3821     labelHashSet entities;
3822
3823     // Handle coupled patches.
3824     handleCoupledPatches(entities);
3825
3826     // Handle layer addition / removal
3827     handleLayerAdditionRemoval();
3828
3829     // Set the thread scheduling sequence
3830     labelList topoSequence(threader_->getNumThreads());
3831
3832     // Linear sequence from 1 to nThreads
3833     forAll(topoSequence, indexI)
3834     {
3835         topoSequence[indexI] = indexI + 1;
3836     }
3837
3838     if (edgeRefinement_)
3839     {
3840         // Initialize stacks
3841         initStacks(entities);
3842
3843         // Execute threads
3844         if (threader_->multiThreaded())
3845         {
3846             executeThreads(topoSequence, handlerPtr_, &edgeRefinementEngine);
3847         }
3848
3849         // Set the master thread to implement modifications
3850         edgeRefinementEngine(&(handlerPtr_[0]));
3851
3852         // Handle mesh slicing events, if necessary
3853         handleMeshSlicing();
3854
3855         if (debug)
3856         {
3857             Info<< nl << "Edge Bisection/Collapse complete." << endl;
3858         }
3859     }
3860
3861     // Re-Initialize stacks
3862     initStacks(entities);
3863
3864     // Execute threads
3865     if (threader_->multiThreaded())
3866     {
3867         if (twoDMesh_)
3868         {
3869             executeThreads(topoSequence, handlerPtr_, &swap2DEdges);
3870         }
3871         else
3872         {
3873             executeThreads(topoSequence, handlerPtr_, &swap3DEdges);
3874         }
3875     }
3876
3877     // Set the master thread to implement modifications
3878     if (twoDMesh_)
3879     {
3880         swap2DEdges(&(handlerPtr_[0]));
3881     }
3882     else
3883     {
3884         swap3DEdges(&(handlerPtr_[0]));
3885     }
3886
3887     if (debug)
3888     {
3889         Info<< nl << "Edge Swapping complete." << endl;
3890     }
3891
3892     // Synchronize coupled patches
3893     syncCoupledPatches(entities);
3894 }
3895
3896
3897 // Reset the mesh and generate mapping information
3898 //  - Return true if topology changes were made.
3899 //  - Return false otherwise (motion only)
3900 bool dynamicTopoFvMesh::resetMesh()
3901 {
3902     // Reduce across processors.
3903     reduce(topoChangeFlag_, orOp<bool>());
3904
3905     if (topoChangeFlag_)
3906     {
3907         // Write out statistics
3908         if (Pstream::parRun() && debug)
3909         {
3910             Pout<< " Bisections :: Total: " << status(3)
3911                 << ", Surface: " << status(5) << nl
3912                 << " Collapses  :: Total: " << status(4)
3913                 << ", Surface: " << status(6) << nl
3914                 << " Swaps      :: Total: " << status(1)
3915                 << ", Surface: " << status(2) << endl;
3916         }
3917         else
3918         {
3919             Info<< " Bisections :: Total: " << status(3)
3920                 << ", Surface: " << status(5) << nl
3921                 << " Collapses  :: Total: " << status(4)
3922                 << ", Surface: " << status(6) << nl
3923                 << " Swaps      :: Total: " << status(1)
3924                 << ", Surface: " << status(2) << endl;
3925         }
3926
3927         if (status(7))
3928         {
3929             Pout<< " Slivers    :: " << status(7) << endl;
3930         }
3931
3932         // Fetch reference to mapper
3933         const topoMapper& fieldMapper = mapper_();
3934
3935         // Set information for the mapping stage
3936         //  - Must be done prior to field-transfers and mesh reset
3937         fieldMapper.storeMeshInformation();
3938
3939         // Set up field-transfers before dealing with mapping
3940         wordList fieldTypes;
3941         List<wordList> fieldNames;
3942         List<List<char> > sendBuffer, recvBuffer;
3943
3944         // Subset fields and transfer
3945         initFieldTransfers
3946         (
3947             fieldTypes,
3948             fieldNames,
3949             sendBuffer,
3950             recvBuffer
3951         );
3952
3953         // Set sizes for mapping
3954         faceWeights_.setSize(facesFromFaces_.size(), scalarField(0));
3955         faceCentres_.setSize(facesFromFaces_.size(), vectorField(0));
3956         cellWeights_.setSize(cellsFromCells_.size(), scalarField(0));
3957         cellCentres_.setSize(cellsFromCells_.size(), vectorField(0));
3958
3959         // Determine if mapping is to be skipped
3960         // Optionally skip mapping for remeshing-only / pre-processing
3961         const dictionary& meshSubDict = dict_.subDict("dynamicTopoFvMesh");
3962
3963         bool skipMapping = false;
3964
3965         if (meshSubDict.found("skipMapping") || mandatory_)
3966         {
3967             skipMapping = readBool(meshSubDict.lookup("skipMapping"));
3968         }
3969
3970         // Fetch the tolerance for mapping
3971         scalar mapTol = 1e-10;
3972
3973         if (meshSubDict.found("mappingTol") || mandatory_)
3974         {
3975             mapTol = readScalar(meshSubDict.lookup("mappingTol"));
3976         }
3977
3978         // Check if outputs are enabled on failure
3979         bool mappingOutput = false;
3980
3981         if (meshSubDict.found("mappingOutput") || mandatory_)
3982         {
3983             mappingOutput = readBool(meshSubDict.lookup("mappingOutput"));
3984         }
3985
3986         clockTime mappingTimer;
3987
3988         // Compute mapping weights for modified entities
3989         threadedMapping(mapTol, skipMapping, mappingOutput);
3990
3991         // Print out stats
3992         Info<< " Mapping time: "
3993             << mappingTimer.elapsedTime() << " s"
3994             << endl;
3995
3996         // Synchronize field transfers prior to the reOrdering stage
3997         syncFieldTransfers
3998         (
3999             fieldTypes,
4000             fieldNames,
4001             recvBuffer
4002         );
4003
4004         // Obtain references to zones, if any
4005         pointZoneMesh& pointZones = polyMesh::pointZones();
4006         faceZoneMesh& faceZones = polyMesh::faceZones();
4007         cellZoneMesh& cellZones = polyMesh::cellZones();
4008
4009         // Allocate temporary lists for mesh-reset
4010         pointField points(nPoints_);
4011         pointField preMotionPoints(nPoints_);
4012         edgeList edges(nEdges_);
4013         faceList faces(nFaces_);
4014         labelList owner(nFaces_);
4015         labelList neighbour(nInternalFaces_);
4016         labelListList faceEdges(nFaces_);
4017         labelListList edgeFaces(nEdges_);
4018         labelList oldPatchStarts(oldPatchStarts_);
4019         labelList oldPatchNMeshPoints(oldPatchNMeshPoints_);
4020         labelListList pointZoneMap(pointZones.size());
4021         labelListList faceZonePointMap(faceZones.size());
4022         labelListList faceZoneFaceMap(faceZones.size());
4023         labelListList cellZoneMap(cellZones.size());
4024
4025         // Obtain faceZone point maps before reordering
4026         List<Map<label> > oldFaceZonePointMaps(faceZones.size());
4027
4028         forAll(faceZones, fzI)
4029         {
4030             oldFaceZonePointMaps[fzI] = faceZones[fzI]().meshPointMap();
4031         }
4032
4033         clockTime reOrderingTimer;
4034
4035         // Reorder the mesh and obtain current topological information
4036         reOrderMesh
4037         (
4038             points,
4039             preMotionPoints,
4040             edges,
4041             faces,
4042             owner,
4043             neighbour,
4044             faceEdges,
4045             edgeFaces,
4046             pointZoneMap,
4047             faceZoneFaceMap,
4048             cellZoneMap
4049         );
4050
4051         // Print out stats
4052         Info<< " Reordering time: "
4053             << reOrderingTimer.elapsedTime() << " s"
4054             << endl;
4055
4056         // Obtain the number of patches before mesh reset
4057         label nOldPatches = boundaryMesh().size();
4058
4059         // Obtain the patch-point maps before resetting the mesh
4060         List<Map<label> > oldPatchPointMaps(nOldPatches);
4061
4062         forAll(oldPatchPointMaps, patchI)
4063         {
4064             oldPatchPointMaps[patchI] = boundaryMesh()[patchI].meshPointMap();
4065         }
4066
4067         // Set weighting information.
4068         // This takes over the weight data.
4069         fieldMapper.setFaceWeights
4070         (
4071             xferMove(faceWeights_),
4072             xferMove(faceCentres_)
4073         );
4074
4075         fieldMapper.setCellWeights
4076         (
4077             xferMove(cellWeights_),
4078             xferMove(cellCentres_)
4079         );
4080
4081         // If the number of patches have changed
4082         // at run-time, reset boundaries first
4083         if (nPatches_ != nOldPatches)
4084         {
4085             resetBoundaries();
4086         }
4087
4088         // Reset the mesh, and specify a non-valid
4089         // boundary to avoid globalData construction
4090         polyMesh::resetPrimitives
4091         (
4092             xferCopy(preMotionPoints),
4093             xferMove(faces),
4094             xferMove(owner),
4095             xferMove(neighbour),
4096             patchSizes_,
4097             patchStarts_,
4098             false
4099         );
4100
4101         // Check the dictionary to determine whether
4102         // edge-connectivity is to be stored on disk.
4103         // This usually benefits restart-time for large
4104         // cases, at the expense of disk-space.
4105         bool storePrimitives = false;
4106
4107         if (meshSubDict.found("storeEdgePrimitives") || mandatory_)
4108         {
4109             storePrimitives =
4110             (
4111                 readBool(meshSubDict.lookup("storeEdgePrimitives"))
4112             );
4113         }
4114
4115         // Reset the edge mesh
4116         eMeshPtr_->resetPrimitives
4117         (
4118             edges,
4119             faceEdges,
4120             edgeFaces,
4121             edgePatchSizes_,
4122             edgePatchStarts_,
4123             true,
4124             (time().outputTime() && storePrimitives)
4125         );
4126
4127         // Generate mapping for points on boundary patches
4128         labelListList patchPointMap(nPatches_);
4129
4130         for (label i = 0; i < nPatches_; i++)
4131         {
4132             // Obtain new patch mesh points after reset.
4133             const labelList& meshPointLabels = boundaryMesh()[i].meshPoints();
4134
4135             patchNMeshPoints_[i] = meshPointLabels.size();
4136
4137             patchPointMap[i].setSize(meshPointLabels.size(), -1);
4138
4139             // Skip for newly introduced patches
4140             if (i >= nOldPatches)
4141             {
4142                 continue;
4143             }
4144
4145             forAll(meshPointLabels, pointI)
4146             {
4147                 label oldIndex = pointMap_[meshPointLabels[pointI]];
4148
4149                 // Check if the point existed before...
4150                 if (oldIndex > -1)
4151                 {
4152                     // Look for the old position on this patch.
4153                     Map<label>::const_iterator oIter =
4154                     (
4155                         oldPatchPointMaps[i].find(oldIndex)
4156                     );
4157
4158                     // Add an entry if the point was found
4159                     if (oIter != oldPatchPointMaps[i].end())
4160                     {
4161                         patchPointMap[i][pointI] = oIter();
4162                     }
4163                 }
4164             }
4165         }
4166
4167         // Generate mapping for faceZone points
4168         forAll(faceZones, fzI)
4169         {
4170             // Obtain new face zone mesh points after reset.
4171             const labelList& meshPointLabels = faceZones[fzI]().meshPoints();
4172
4173             faceZonePointMap[fzI].setSize(meshPointLabels.size(), -1);
4174
4175             forAll(meshPointLabels, pointI)
4176             {
4177                 label oldIndex = pointMap_[meshPointLabels[pointI]];
4178
4179                 // Check if the point existed before...
4180                 if (oldIndex > -1)
4181                 {
4182                     // Look for the old position on this patch.
4183                     Map<label>::const_iterator oIter =
4184                     (
4185                         oldFaceZonePointMaps[fzI].find(oldIndex)
4186                     );
4187
4188                     // Add an entry if the point was found
4189                     if (oIter != oldFaceZonePointMaps[fzI].end())
4190                     {
4191                         faceZonePointMap[fzI][pointI] = oIter();
4192                     }
4193                 }
4194             }
4195         }
4196
4197         // Generate new mesh mapping information
4198         mapPolyMesh mpm
4199         (
4200             (*this),
4201             nOldPoints_,
4202             nOldFaces_,
4203             nOldCells_,
4204             pointMap_,
4205             pointsFromPoints_,
4206             faceMap_,
4207             facesFromPoints_,
4208             facesFromEdges_,
4209             facesFromFaces_,
4210             cellMap_,
4211             cellsFromPoints_,
4212             cellsFromEdges_,
4213             cellsFromFaces_,
4214             cellsFromCells_,
4215             reversePointMap_,
4216             reverseFaceMap_,
4217             reverseCellMap_,
4218             flipFaces_,
4219             patchPointMap,
4220             pointZoneMap,
4221             faceZonePointMap,
4222             faceZoneFaceMap,
4223             cellZoneMap,
4224             preMotionPoints,
4225             oldPatchStarts,
4226             oldPatchNMeshPoints,
4227             true
4228         );
4229
4230         // Update the underlying mesh, and map all related fields
4231         updateMesh(mpm);
4232
4233         if (mpm.hasMotionPoints())
4234         {
4235             // Perform a dummy movePoints to force V0 creation
4236             movePoints(mpm.preMotionPoints());
4237
4238             // Reset old-volumes
4239             resetMotion();
4240             setV0();
4241         }
4242
4243         // Correct volume fluxes on the old mesh
4244         fieldMapper.correctFluxes();
4245
4246         // Clear mapper after use
4247         fieldMapper.clear();
4248
4249         if (mpm.hasMotionPoints())
4250         {
4251             // Now move mesh to new points and
4252             // compute correct mesh-fluxes.
4253             movePoints(points);
4254         }
4255
4256         // Update the mesh-motion solver
4257         if (motionSolver_.valid())
4258         {
4259             motionSolver_->updateMesh(mpm);
4260         }
4261
4262         // Clear unwanted member data
4263         addedFacePatches_.clear();
4264         addedEdgePatches_.clear();
4265         addedPointZones_.clear();
4266         addedFaceZones_.clear();
4267         addedCellZones_.clear();
4268         faceParents_.clear();
4269         cellParents_.clear();
4270
4271         // Clear the deleted entity map
4272         deletedPoints_.clear();
4273         deletedEdges_.clear();
4274         deletedFaces_.clear();
4275         deletedCells_.clear();
4276
4277         // Clear flipFaces
4278         flipFaces_.clear();
4279
4280         // Set new sizes for the reverse maps
4281         reversePointMap_.setSize(nPoints_, -7);
4282         reverseEdgeMap_.setSize(nEdges_, -7);
4283         reverseFaceMap_.setSize(nFaces_, -7);
4284         reverseCellMap_.setSize(nCells_, -7);
4285
4286         // Update "old" information
4287         nOldPoints_ = nPoints_;
4288         nOldEdges_ = nEdges_;
4289         nOldFaces_ = nFaces_;
4290         nOldCells_ = nCells_;
4291         nOldInternalFaces_ = nInternalFaces_;
4292         nOldInternalEdges_ = nInternalEdges_;
4293
4294         for (label i = 0; i < nPatches_; i++)
4295         {
4296             oldPatchSizes_[i] = patchSizes_[i];
4297             oldPatchStarts_[i] = patchStarts_[i];
4298             oldEdgePatchSizes_[i] = edgePatchSizes_[i];
4299             oldEdgePatchStarts_[i] = edgePatchStarts_[i];
4300             oldPatchNMeshPoints_[i] = patchNMeshPoints_[i];
4301         }
4302
4303         // Clear parallel structures
4304         if (Pstream::parRun())
4305         {
4306             procIndices_.clear();
4307             sendMeshes_.clear();
4308             recvMeshes_.clear();
4309         }
4310
4311         bool checkCplBoundaries = false;
4312
4313         if (meshSubDict.found("checkCoupledBoundaries") || mandatory_)
4314         {
4315             checkCplBoundaries =
4316             (
4317                 readBool(meshSubDict.lookup("checkCoupledBoundaries"))
4318             );
4319         }
4320
4321         if (checkCplBoundaries)
4322         {
4323             bool failed = checkCoupledBoundaries();
4324
4325             if (failed)
4326             {
4327                 FatalErrorIn("bool dynamicTopoFvMesh::resetMesh()")
4328                     << " Coupled boundary check failed on processor: "
4329                     << Pstream::myProcNo()
4330                     << abort(FatalError);
4331             }
4332         }
4333
4334         // Reset statistics
4335         statistics_ = 0;
4336     }
4337     else
4338     {
4339         // No topology changes were made.
4340         // Only execute mesh-motion.
4341         if (motionSolver_.valid())
4342         {
4343             movePoints(points_);
4344         }
4345     }
4346
4347     // Obtain mesh stats after topo-changes
4348     meshQuality(true);
4349
4350     // Dump length-scale to disk, if requested.
4351     calculateLengthScale(true);
4352
4353     // Reset and return flag
4354     if (topoChangeFlag_)
4355     {
4356         topoChangeFlag_ = false;
4357         return true;
4358     }
4359
4360     // No changes were made.
4361     return false;
4362 }
4363
4364
4365 // Update mesh corresponding to the given map
4366 void dynamicTopoFvMesh::updateMesh(const mapPolyMesh& mpm)
4367 {
4368     if (coupledModification_)
4369     {
4370         // This bit gets called only during the load-balancing
4371         // stage, since the fvMesh::updateMesh is a bit different
4372         fvMesh::updateMesh(mpm);
4373         return;
4374     }
4375
4376     // Delete oldPoints in polyMesh
4377     polyMesh::resetMotion();
4378
4379     // Clear-out fvMesh geometry and addressing
4380     fvMesh::clearOut();
4381
4382     // Update polyMesh.
4383     polyMesh::updateMesh(mpm);
4384
4385     // Map all fields
4386     mapFields(mpm);
4387 }
4388
4389
4390 // Map all fields in time using a customized mapper
4391 void dynamicTopoFvMesh::mapFields(const mapPolyMesh& mpm) const
4392 {
4393     if (coupledModification_)
4394     {
4395         // This bit gets called only during the load-balancing
4396         // stage, since the fvMesh::mapFields is a bit different
4397         fvMesh::mapFields(mpm);
4398         return;
4399     }
4400
4401     if (debug)
4402     {
4403         Info<< "void dynamicTopoFvMesh::mapFields(const mapPolyMesh&) const: "
4404             << "Mapping fv fields."
4405             << endl;
4406     }
4407
4408     const topoMapper& fieldMapper = mapper_();
4409
4410     // Set the mapPolyMesh object in the mapper
4411     fieldMapper.setMapper(mpm);
4412
4413     // Conservatively map scalar/vector volFields
4414     conservativeMapVolFields<scalar>(fieldMapper);
4415     conservativeMapVolFields<vector>(fieldMapper);
4416
4417     // Map all the volFields in the objectRegistry
4418     MapGeometricFields<sphericalTensor,fvPatchField,topoMapper,volMesh>
4419         (fieldMapper);
4420     MapGeometricFields<symmTensor,fvPatchField,topoMapper,volMesh>
4421         (fieldMapper);
4422     MapGeometricFields<tensor,fvPatchField,topoMapper,volMesh>
4423         (fieldMapper);
4424
4425     // Conservatively map scalar/vector surfaceFields
4426     conservativeMapSurfaceFields<scalar>(fieldMapper);
4427     conservativeMapSurfaceFields<vector>(fieldMapper);
4428
4429     // Map all the surfaceFields in the objectRegistry
4430     MapGeometricFields<sphericalTensor,fvsPatchField,topoMapper,surfaceMesh>
4431         (fieldMapper);
4432     MapGeometricFields<symmTensor,fvsPatchField,topoMapper,surfaceMesh>
4433         (fieldMapper);
4434     MapGeometricFields<tensor,fvsPatchField,topoMapper,surfaceMesh>
4435         (fieldMapper);
4436 }
4437
4438
4439 // Update the mesh for motion / topology changes.
4440 //  - Return true if topology changes have occurred
4441 bool dynamicTopoFvMesh::update()
4442 {
4443     // Re-read options, in case they have been modified at run-time
4444     readOptionalParameters(true);
4445
4446     // Set old point positions
4447     oldPoints_ = polyMesh::points();
4448
4449     // Invoke mesh-motion solver and store new points
4450     if (motionSolver_.valid())
4451     {
4452         points_ = motionSolver_->newPoints()();
4453     }
4454     else
4455     {
4456         // Set point positions from mesh
4457         points_ = polyMesh::points();
4458     }
4459
4460     // Obtain mesh stats before topo-changes
4461     bool noSlivers = meshQuality(true);
4462
4463     // Return if the interval is invalid,
4464     // not at re-mesh interval, or slivers are absent.
4465     // Handy while using only mesh-motion.
4466     if (interval_ < 0 || ((time().timeIndex() % interval_ != 0) && noSlivers))
4467     {
4468         return resetMesh();
4469     }
4470
4471     // Calculate the edge length-scale for the mesh
4472     calculateLengthScale();
4473
4474     // Track mesh topology modification time
4475     clockTime topoTimer;
4476
4477     // Invoke the threaded topoModifier
4478     threadedTopoModifier();
4479
4480     Info<< " Topo modifier time: "
4481         << topoTimer.elapsedTime() << " s"
4482         << endl;
4483
4484     // Apply all topology changes (if any) and reset mesh.
4485     return resetMesh();
4486 }
4487
4488
4489 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
4490
4491 void dynamicTopoFvMesh::operator=(const dynamicTopoFvMesh& rhs)
4492 {
4493     // Check for assignment to self
4494     if (this == &rhs)
4495     {
4496         FatalErrorIn
4497         (
4498             "dynamicTopoFvMesh::operator=(const dynamicTopoFvMesh&)"
4499         )
4500             << "Attempted assignment to self"
4501             << abort(FatalError);
4502     }
4503 }
4504
4505
4506 } // End namespace Foam
4507
4508 // ************************************************************************* //