src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/dynamicTopoFvMeshCheck.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     Functions specific to connectivity checking and debugging
  30
  31 Author
  32     Sandeep Menon
  33     University of Massachusetts Amherst
  34     All rights reserved
  35
  36 \*---------------------------------------------------------------------------*/
  37
  38 #include "dynamicTopoFvMesh.H"
  39
  40 #include "IOmanip.H"
  41 #include "volFields.H"
  42 #include "triPointRef.H"
  43 #include "tetPointRef.H"
  44 #include "coupledInfo.H"
  45
  46 namespace Foam
  47 {
  48
  49 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  50
  51 // Compute mesh-quality, and return true if no slivers are present
  52 bool dynamicTopoFvMesh::meshQuality
  53 (
  54     bool outputOption
  55 )
  56 {
  57     // Compute statistics on the fly
  58     label nCells = 0, minCell = -1;
  59     scalar maxQuality = -GREAT;
  60     scalar minQuality =  GREAT;
  61     scalar cQuality = 0.0;
  62     scalar meanQuality = 0.0;
  63
  64     // Track slivers
  65     bool sliversAbsent = true;
  66     thresholdSlivers_.clear();
  67
  68     // Loop through all cells in the mesh and compute cell quality
  69     forAll(cells_, cellI)
  70     {
  71         const cell& cellToCheck = cells_[cellI];
  72
  73         if (cellToCheck.empty())
  74         {
  75             continue;
  76         }
  77
  78         // Skip hexahedral cells
  79         if (cellToCheck.size() == 6)
  80         {
  81             cQuality = 1.0;
  82             meanQuality += cQuality;
  83
  84             // Update min / max
  85             maxQuality = Foam::max(cQuality, maxQuality);
  86             minQuality = Foam::min(cQuality, minQuality);
  87
  88             nCells++;
  89
  90             continue;
  91         }
  92
  93         if (twoDMesh_)
  94         {
  95             // Assume XY plane here
  96             vector n = vector(0,0,1);
  97
  98             // Get a triangular boundary face
  99             forAll(cellToCheck, faceI)
 100             {
 101                 const face& faceToCheck = faces_[cellToCheck[faceI]];
 102
 103                 if (faceToCheck.size() == 3)
 104                 {
 105                     triPointRef tpr
 106                     (
 107                         points_[faceToCheck[0]],
 108                         points_[faceToCheck[1]],
 109                         points_[faceToCheck[2]]
 110                     );
 111
 112                     // Assume centre-plane passes through origin
 113                     cQuality =
 114                     (
 115                         tpr.quality() *
 116                         (
 117                             Foam::sign
 118                             (
 119                                 tpr.normal() &
 120                                 ((tpr.centre() & n) * n)
 121                             )
 122                         )
 123                     );
 124
 125                     break;
 126                 }
 127             }
 128         }
 129         else
 130         {
 131             const label bfIndex = cellToCheck[0];
 132             const label cfIndex = cellToCheck[1];
 133
 134             const face& baseFace = faces_[bfIndex];
 135             const face& checkFace = faces_[cfIndex];
 136
 137             // Get the fourth point
 138             label apexPoint =
 139             (
 140                 meshOps::findIsolatedPoint(baseFace, checkFace)
 141             );
 142
 143             // Compute cell quality
 144             if (owner_[bfIndex] == cellI)
 145             {
 146                 cQuality =
 147                 (
 148                     tetMetric_
 149                     (
 150                         points_[baseFace[2]],
 151                         points_[baseFace[1]],
 152                         points_[baseFace[0]],
 153                         points_[apexPoint]
 154                     )
 155                 );
 156             }
 157             else
 158             {
 159                 cQuality =
 160                 (
 161                     tetMetric_
 162                     (
 163                         points_[baseFace[0]],
 164                         points_[baseFace[1]],
 165                         points_[baseFace[2]],
 166                         points_[apexPoint]
 167                     )
 168                 );
 169             }
 170         }
 171
 172         // Update statistics
 173         maxQuality = Foam::max(cQuality, maxQuality);
 174
 175         if (cQuality < minQuality)
 176         {
 177             minQuality = cQuality;
 178             minCell = cellI;
 179         }
 180
 181         meanQuality += cQuality;
 182         nCells++;
 183
 184         // Add to the list of slivers
 185         if ((cQuality < sliverThreshold_) && (cQuality > 0.0))
 186         {
 187             thresholdSlivers_.insert(cellI, cQuality);
 188         }
 189     }
 190
 191     if (thresholdSlivers_.size())
 192     {
 193         sliversAbsent = false;
 194     }
 195
 196     // Reduce across processors.
 197     reduce(sliversAbsent, andOp<bool>());
 198
 199     // Output statistics:
 200     if (outputOption || (debug > 0))
 201     {
 202         // Reduce statistics across processors.
 203         reduce(minQuality, minOp<scalar>());
 204         reduce(maxQuality, maxOp<scalar>());
 205         reduce(meanQuality, sumOp<scalar>());
 206         reduce(nCells, sumOp<label>());
 207
 208         if (minQuality < 0.0)
 209         {
 210             WarningIn
 211             (
 212                 "bool dynamicTopoFvMesh::meshQuality"
 213                 "(bool outputOption)"
 214             )
 215                 << nl
 216                 << "Minimum cell quality is: " << minQuality
 217                 << " at cell: " << minCell
 218                 << endl;
 219         }
 220
 221         Info<< nl
 222             << " ~~~ Mesh Quality Statistics ~~~ " << nl
 223             << " Min: " << minQuality << nl
 224             << " Max: " << maxQuality << nl
 225             << " Mean: " << meanQuality/nCells << nl
 226             << " Cells: " << nCells << nl
 227             << " ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ " << nl
 228             << endl;
 229     }
 230
 231     return sliversAbsent;
 232 }
 233
 234
 235 // Utility to check whether points of an edge lie on a boundary.
 236 const FixedList<bool,2>
 237 dynamicTopoFvMesh::checkEdgeBoundary
 238 (
 239     const label eIndex
 240 ) const
 241 {
 242     FixedList<bool,2> edgeBoundary(false);
 243
 244     const edge& edgeToCheck = edges_[eIndex];
 245
 246     // Loop through edges connected to both points,
 247     // and check if any of them lie on boundaries.
 248     // Used to ensure that collapses happen towards boundaries.
 249     forAll(edgeToCheck, pointI)
 250     {
 251         const labelList& pEdges = pointEdges_[edgeToCheck[pointI]];
 252
 253         forAll(pEdges, edgeI)
 254         {
 255             // Determine the patch this edge belongs to
 256             if (whichEdgePatch(pEdges[edgeI]) > -1)
 257             {
 258                 edgeBoundary[pointI] = true;
 259                 break;
 260             }
 261         }
 262     }
 263
 264     return edgeBoundary;
 265 }
 266
 267
 268 // Check whether the given edge is on a bounding curve
 269 //  - If nProcCurves is provided, the variable is incremented
 270 //    if the edge is processor-coupled
 271 bool dynamicTopoFvMesh::checkBoundingCurve
 272 (
 273     const label eIndex,
 274     const bool overRidePurityCheck,
 275     label* nProcCurves
 276 ) const
 277 {
 278     // Internal edges don't count
 279     label edgePatch = -1;
 280
 281     // If this entity was deleted, skip it.
 282     if (edgeFaces_[eIndex].empty())
 283     {
 284         // Return true so that swap3DEdges skips this edge.
 285         return true;
 286     }
 287
 288     // Check if two boundary faces lie on different face-patches
 289     bool procCoupled = false;
 290     FixedList<label, 2> fPatches(-1);
 291     FixedList<vector, 2> fNorm(vector::zero);
 292
 293     if ((edgePatch = whichEdgePatch(eIndex)) < 0)
 294     {
 295         return false;
 296     }
 297     else
 298     {
 299         // Check whether this edge shouldn't be swapped
 300         if (findIndex(noSwapPatchIDs_, edgePatch) > -1)
 301         {
 302             return true;
 303         }
 304
 305         // Explicit check for processor edges (both 2D and 3D)
 306         if (processorCoupledEntity(eIndex, false, true))
 307         {
 308             // Increment nProcCurves
 309             if (nProcCurves)
 310             {
 311                 (*nProcCurves)++;
 312             }
 313
 314             // Check for pure processor edge, and if not,
 315             // fetch boundary patch labels / normals
 316             if
 317             (
 318                 processorCoupledEntity
 319                 (
 320                     eIndex,
 321                     false,
 322                     true,
 323                     true,
 324                     &fPatches,
 325                     &fNorm
 326                 )
 327             )
 328             {
 329                 // 'Pure' processor coupled edges don't count
 330                 return false;
 331             }
 332             else
 333             if (!overRidePurityCheck)
 334             {
 335                 // This edge lies between a processor and physical patch,
 336                 //  - This a bounding curve (unless an override is requested)
 337                 //  - An override is warranted for 2-2 swaps on impure edges,
 338                 //    which is typically requested by swap3DEdges.
 339                 return true;
 340             }
 341
 342             // Specify that the edge is procCoupled
 343             procCoupled = true;
 344         }
 345     }
 346
 347     if (procCoupled)
 348     {
 349         // Normalize patch normals from coupled check
 350         fNorm[0] /= mag(fNorm[0]) + VSMALL;
 351         fNorm[1] /= mag(fNorm[1]) + VSMALL;
 352     }
 353     else
 354     {
 355         // Fetch patch indices / normals
 356         const labelList& eFaces = edgeFaces_[eIndex];
 357
 358         label fPatch = -1, count = 0;
 359
 360         forAll(eFaces, faceI)
 361         {
 362             if ((fPatch = whichPatch(eFaces[faceI])) > -1)
 363             {
 364                 // Obtain the normal.
 365                 fNorm[count] = faces_[eFaces[faceI]].normal(points_);
 366
 367                 // Normalize it.
 368                 fNorm[count] /= mag(fNorm[count]) + VSMALL;
 369
 370                 // Note patch index
 371                 fPatches[count] = fPatch;
 372
 373                 count++;
 374
 375                 if (count == 2)
 376                 {
 377                     break;
 378                 }
 379             }
 380         }
 381     }
 382
 383     // Check for legitimate patches
 384     if (fPatches[0] < 0 || fPatches[1] < 0)
 385     {
 386         const labelList& eFaces = edgeFaces_[eIndex];
 387
 388         forAll(eFaces, faceI)
 389         {
 390             Pout<< " Face: " << eFaces[faceI]
 391                 << " :: " << faces_[eFaces[faceI]]
 392                 << " Patch: " << whichPatch(eFaces[faceI]) << nl;
 393         }
 394
 395         label epI = whichEdgePatch(eIndex);
 396
 397         FatalErrorIn
 398         (
 399             "bool dynamicTopoFvMesh::checkBoundingCurve"
 400             "(const label, const bool) const"
 401         )
 402             << " Edge: " << eIndex << ":: " << edges_[eIndex]
 403             << " Patch: "
 404             << (epI < 0 ? "Internal" : boundaryMesh()[epI].name())
 405             << " edgeFaces: " << eFaces << nl
 406             << " expected 2 boundary patches." << nl
 407             << " fPatches[0]: " << fPatches[0] << nl
 408             << " fPatches[1]: " << fPatches[1] << nl
 409             << " fNorm[0]: " << fNorm[0] << nl
 410             << " fNorm[1]: " << fNorm[1] << nl
 411             << " coupledModification: " << coupledModification_
 412             << abort(FatalError);
 413     }
 414
 415     scalar deviation = (fNorm[0] & fNorm[1]);
 416
 417     // Check if the swap-curvature is too high
 418     if (mag(deviation) < swapDeviation_)
 419     {
 420         return true;
 421     }
 422
 423     // Check if the edge borders two different patches
 424     if (fPatches[0] != fPatches[1])
 425     {
 426         return true;
 427     }
 428
 429     // Not on a bounding curve
 430     return false;
 431 }
 432
 433
 434 // Check triangulation quality for an edge index
 435 bool dynamicTopoFvMesh::checkQuality
 436 (
 437     const label eIndex,
 438     const labelList& m,
 439     const PtrList<scalarListList>& Q,
 440     const scalar minQuality,
 441     const label checkIndex
 442 ) const
 443 {
 444     bool myResult = false;
 445
 446     // Non-coupled check
 447     if (Q[checkIndex][0][m[checkIndex]-1] > minQuality)
 448     {
 449         myResult = true;
 450
 451         if (debug > 2)
 452         {
 453             Pout<< nl << nl
 454                 << " eIndex: " << eIndex
 455                 << " minQuality: " << minQuality
 456                 << " newQuality: " << Q[checkIndex][0][m[checkIndex]-1]
 457                 << endl;
 458         }
 459     }
 460
 461     if (coupledModification_)
 462     {
 463         // Only locally coupled indices require checks
 464         if (locallyCoupledEntity(eIndex))
 465         {
 466             // Check the quality of the slave edge as well.
 467             label sIndex = -1;
 468
 469             // Loop through masterToSlave and determine the slave index.
 470             forAll(patchCoupling_, patchI)
 471             {
 472                 if (patchCoupling_(patchI))
 473                 {
 474                     const label edgeEnum  = coupleMap::EDGE;
 475                     const coupleMap& cMap = patchCoupling_[patchI].map();
 476
 477                     if ((sIndex = cMap.findSlave(edgeEnum, eIndex)) > -1)
 478                     {
 479                         break;
 480                     }
 481                 }
 482             }
 483
 484             if (sIndex == -1)
 485             {
 486                 FatalErrorIn
 487                 (
 488                     "bool dynamicTopoFvMesh::checkQuality\n"
 489                     "(\n"
 490                     "    const label eIndex,\n"
 491                     "    const labelList& m,\n"
 492                     "    const PtrList<scalarListList>& Q,\n"
 493                     "    const scalar minQuality,\n"
 494                     "    const label checkIndex\n"
 495                     ") const\n"
 496                 )
 497                     << "Coupled maps were improperly specified." << nl
 498                     << " Slave index not found for: " << nl
 499                     << " Edge: " << eIndex << nl
 500                     << abort(FatalError);
 501             }
 502
 503             // Turn off switch temporarily.
 504             unsetCoupledModification();
 505
 506             // Recursively call for the slave edge.
 507             myResult =
 508             (
 509                 myResult && checkQuality(sIndex, m, Q, minQuality, 1)
 510             );
 511
 512             // Turn it back on.
 513             setCoupledModification();
 514         }
 515     }
 516
 517     return myResult;
 518 }
 519
 520
 521 // Print out tables for debugging
 522 void dynamicTopoFvMesh::printTables
 523 (
 524     const labelList& m,
 525     const PtrList<scalarListList>& Q,
 526     const PtrList<labelListList>& K,
 527     const label checkIndex
 528 ) const
 529 {
 530     Pout<< "m: " << m[checkIndex] << endl;
 531
 532     // Print out Q
 533     Pout<< "===" << nl
 534         << " Q " << nl
 535         << "===" << endl;
 536
 537     Pout<< "   ";
 538
 539     for (label j = 0; j < m[checkIndex]; j++)
 540     {
 541         Pout<< setw(12) << j;
 542     }
 543
 544     Pout<< nl;
 545
 546     for (label j = 0; j < m[checkIndex]; j++)
 547     {
 548         Pout<< "-------------";
 549     }
 550
 551     Pout<< nl;
 552
 553     for (label i = 0; i < (m[checkIndex]-2); i++)
 554     {
 555         Pout<< i << ": ";
 556
 557         for (label j = 0; j < m[checkIndex]; j++)
 558         {
 559             Pout<< setw(12) << Q[checkIndex][i][j];
 560         }
 561
 562         Pout<< nl;
 563     }
 564
 565     // Print out K
 566     Pout<< "===" << nl
 567         << " K " << nl
 568         << "===" << endl;
 569
 570     Pout<< "   ";
 571
 572     for (label j = 0; j < m[checkIndex]; j++)
 573     {
 574         Pout<< setw(12) << j;
 575     }
 576
 577     Pout<< nl;
 578
 579     for (label i = 0; i < (m[checkIndex]-2); i++)
 580     {
 581         Pout<< i << ": ";
 582
 583         for (label j = 0; j < m[checkIndex]; j++)
 584         {
 585             Pout<< setw(12) << K[checkIndex][i][j];
 586         }
 587
 588         Pout<< nl;
 589     }
 590
 591     Pout<< endl;
 592 }
 593
 594
 595 // Check old-volumes for an input triangulation
 596 bool dynamicTopoFvMesh::checkTriangulationVolumes
 597 (
 598     const label eIndex,
 599     const labelList& hullVertices,
 600     const labelListList& triangulations
 601 ) const
 602 {
 603     label m = hullVertices.size();
 604     scalar oldTetVol = 0.0;
 605     scalar newTetVol = 0.0;
 606
 607     const edge& edgeToCheck = edges_[eIndex];
 608
 609     for (label i = 0; i < (m-2); i++)
 610     {
 611         // Compute volume for the upper-half
 612         newTetVol =
 613         (
 614             tetPointRef
 615             (
 616                 points_[hullVertices[triangulations[0][i]]],
 617                 points_[hullVertices[triangulations[1][i]]],
 618                 points_[hullVertices[triangulations[2][i]]],
 619                 points_[edgeToCheck[0]]
 620             ).mag()
 621         );
 622
 623         oldTetVol =
 624         (
 625             tetPointRef
 626             (
 627                 oldPoints_[hullVertices[triangulations[0][i]]],
 628                 oldPoints_[hullVertices[triangulations[1][i]]],
 629                 oldPoints_[hullVertices[triangulations[2][i]]],
 630                 oldPoints_[edgeToCheck[0]]
 631             ).mag()
 632         );
 633
 634         if (oldTetVol < 0.0 || (mag(oldTetVol) < mag(0.1 * newTetVol)))
 635         {
 636             if (debug > 2)
 637             {
 638                 Pout<< " Swap sequence leads to bad old-volumes." << nl
 639                     << " Edge: " << edgeToCheck << nl
 640                     << " using Points: " << nl
 641                     << oldPoints_[hullVertices[triangulations[0][i]]] << nl
 642                     << oldPoints_[hullVertices[triangulations[1][i]]] << nl
 643                     << oldPoints_[hullVertices[triangulations[2][i]]] << nl
 644                     << oldPoints_[edgeToCheck[0]] << nl
 645                     << " Old Volume: " << oldTetVol << nl
 646                     << " New Volume: " << newTetVol << nl
 647                     << endl;
 648             }
 649
 650             return true;
 651         }
 652
 653         newTetVol =
 654         (
 655             tetPointRef
 656             (
 657                 points_[hullVertices[triangulations[2][i]]],
 658                 points_[hullVertices[triangulations[1][i]]],
 659                 points_[hullVertices[triangulations[0][i]]],
 660                 points_[edgeToCheck[1]]
 661             ).mag()
 662         );
 663
 664         oldTetVol =
 665         (
 666             tetPointRef
 667             (
 668                 oldPoints_[hullVertices[triangulations[2][i]]],
 669                 oldPoints_[hullVertices[triangulations[1][i]]],
 670                 oldPoints_[hullVertices[triangulations[0][i]]],
 671                 oldPoints_[edgeToCheck[1]]
 672             ).mag()
 673         );
 674
 675         if (oldTetVol < 0.0 || (mag(oldTetVol) < mag(0.1 * newTetVol)))
 676         {
 677             if (debug > 2)
 678             {
 679                 Pout<< " Swap sequence leads to bad old-volumes." << nl
 680                     << " Edge: " << edgeToCheck << nl
 681                     << " using Points: " << nl
 682                     << oldPoints_[hullVertices[triangulations[2][i]]] << nl
 683                     << oldPoints_[hullVertices[triangulations[1][i]]] << nl
 684                     << oldPoints_[hullVertices[triangulations[0][i]]] << nl
 685                     << oldPoints_[edgeToCheck[1]] << nl
 686                     << " Old Volume: " << oldTetVol << nl
 687                     << " New Volume: " << newTetVol << nl
 688                     << endl;
 689             }
 690
 691             return true;
 692         }
 693     }
 694
 695     return false;
 696 }
 697
 698
 699 // Write out connectivity for an edge
 700 void dynamicTopoFvMesh::writeEdgeConnectivity
 701 (
 702     const label eIndex
 703 ) const
 704 {
 705     // Write out edge
 706     writeVTK("Edge_" + Foam::name(eIndex), eIndex, 1, false, true);
 707
 708     const labelList& eFaces = edgeFaces_[eIndex];
 709
 710     // Write out edge faces
 711     writeVTK
 712     (
 713         "EdgeFaces_"
 714       + Foam::name(eIndex)
 715       + '_'
 716       + Foam::name(Pstream::myProcNo()),
 717         eFaces,
 718         2, false, true
 719     );
 720
 721     DynamicList<label> edgeCells(10);
 722
 723     forAll(eFaces, faceI)
 724     {
 725         label pIdx = whichPatch(eFaces[faceI]);
 726
 727         word pName((pIdx < 0) ? "Internal" : boundaryMesh()[pIdx].name());
 728
 729         Pout<< " Face: " << eFaces[faceI]
 730             << " :: " << faces_[eFaces[faceI]]
 731             << " Patch: " << pName
 732             << " Proc: " << Pstream::myProcNo() << nl;
 733
 734         label own = owner_[eFaces[faceI]];
 735         label nei = neighbour_[eFaces[faceI]];
 736
 737         if (findIndex(edgeCells, own) == -1)
 738         {
 739             edgeCells.append(own);
 740         }
 741
 742         if (nei == -1)
 743         {
 744             continue;
 745         }
 746
 747         if (findIndex(edgeCells, nei) == -1)
 748         {
 749             edgeCells.append(nei);
 750         }
 751     }
 752
 753     // Write out cells connected to edge
 754     writeVTK
 755     (
 756         "EdgeCells_"
 757       + Foam::name(eIndex)
 758       + '_'
 759       + Foam::name(Pstream::myProcNo()),
 760         edgeCells,
 761         3, false, true
 762     );
 763
 764     if (twoDMesh_)
 765     {
 766         return;
 767     }
 768
 769     // Check processors for coupling
 770     forAll(procIndices_, pI)
 771     {
 772         // Fetch reference to subMesh
 773         const coupleMap& cMap = recvMeshes_[pI].map();
 774         const dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
 775
 776         label sI = -1;
 777
 778         if ((sI = cMap.findSlave(coupleMap::EDGE, eIndex)) == -1)
 779         {
 780             continue;
 781         }
 782
 783         const edge& slaveEdge = mesh.edges_[sI];
 784         const labelList& seFaces = mesh.edgeFaces_[sI];
 785
 786         edge cE
 787         (
 788             cMap.findMaster(coupleMap::POINT, slaveEdge[0]),
 789             cMap.findMaster(coupleMap::POINT, slaveEdge[1])
 790         );
 791
 792         Pout<< " >> Edge: " << sI << "::" << slaveEdge
 793             << " mapped: " << cE << nl;
 794
 795         mesh.writeVTK
 796         (
 797             "EdgeFaces_"
 798           + Foam::name(eIndex)
 799           + '_'
 800           + Foam::name(procIndices_[pI]),
 801             seFaces,
 802             2, false, true
 803         );
 804
 805         // Clear existing list
 806         edgeCells.clear();
 807
 808         forAll(seFaces, faceI)
 809         {
 810             label pIdx = mesh.whichPatch(seFaces[faceI]);
 811
 812             word pName
 813             (
 814                 (pIdx < 0) ?
 815                 "Internal" :
 816                 mesh.boundaryMesh()[pIdx].name()
 817             );
 818
 819             Pout<< " Face: " << seFaces[faceI]
 820                 << " :: " << mesh.faces_[seFaces[faceI]]
 821                 << " Patch: " << pName
 822                 << " Proc: " << procIndices_[pI] << nl;
 823
 824             label own = mesh.owner_[seFaces[faceI]];
 825             label nei = mesh.neighbour_[seFaces[faceI]];
 826
 827             if (findIndex(edgeCells, own) == -1)
 828             {
 829                 edgeCells.append(own);
 830             }
 831
 832             if (nei == -1)
 833             {
 834                 continue;
 835             }
 836
 837             if (findIndex(edgeCells, nei) == -1)
 838             {
 839                 edgeCells.append(nei);
 840             }
 841         }
 842
 843         mesh.writeVTK
 844         (
 845             "EdgeCells_"
 846           + Foam::name(eIndex)
 847           + '_'
 848           + Foam::name(procIndices_[pI]),
 849             edgeCells,
 850             3, false, true
 851         );
 852     }
 853 }
 854
 855
 856 // Output an entity as a VTK file
 857 void dynamicTopoFvMesh::writeVTK
 858 (
 859     const word& name,
 860     const label entity,
 861     const label primitiveType,
 862     const bool useOldConnectivity,
 863     const bool useOldPoints
 864 ) const
 865 {
 866     writeVTK
 867     (
 868         name,
 869         labelList(1, entity),
 870         primitiveType,
 871         useOldConnectivity,
 872         useOldPoints
 873     );
 874 }
 875
 876
 877 // Output a list of primitives as a VTK file.
 878 //  - primitiveType is:
 879 //      0: List of points
 880 //      1: List of edges
 881 //      2: List of faces
 882 //      3: List of cells
 883 void dynamicTopoFvMesh::writeVTK
 884 (
 885     const word& name,
 886     const labelList& cList,
 887     const label primitiveType,
 888     const bool useOldConnectivity,
 889     const bool useOldPoints,
 890     const UList<scalar>& scalField,
 891     const UList<label>& lablField
 892 ) const
 893 {
 894     // Check if spatial bounding box has been specified
 895     const dictionary& meshSubDict = dict_.subDict("dynamicTopoFvMesh");
 896
 897     labelList entityList;
 898
 899     if (meshSubDict.found("spatialDebug") && !useOldConnectivity)
 900     {
 901         // Read the bounding box
 902         boundBox bb
 903         (
 904             meshSubDict.subDict("spatialDebug").lookup("debugBoundBox")
 905         );
 906
 907         DynamicList<label> cSubList(10);
 908
 909         forAll(cList, cellI)
 910         {
 911             label index = cList[cellI];
 912
 913             if (index < 0)
 914             {
 915                 continue;
 916             }
 917
 918             point containPoint(vector::zero);
 919
 920             switch (primitiveType)
 921             {
 922                 // Are we looking at points?
 923                 case 0:
 924                 {
 925                     containPoint = points_[index];
 926                     break;
 927                 }
 928
 929                 // Are we looking at edges?
 930                 case 1:
 931                 {
 932                     containPoint = edges_[index].centre(points_);
 933                     break;
 934                 }
 935
 936                 // Are we looking at faces?
 937                 case 2:
 938                 {
 939                     containPoint = faces_[index].centre(points_);
 940                     break;
 941                 }
 942
 943                 // Are we looking at cells?
 944                 case 3:
 945                 {
 946                     scalar volume = 0.0;
 947
 948                     // Compute centre
 949                     meshOps::cellCentreAndVolume
 950                     (
 951                         index,
 952                         points_,
 953                         faces_,
 954                         cells_,
 955                         owner_,
 956                         containPoint,
 957                         volume
 958                     );
 959
 960                     break;
 961                 }
 962             }
 963
 964             // Is the point of interest?
 965             if (bb.contains(containPoint))
 966             {
 967                 cSubList.append(index);
 968             }
 969         }
 970
 971         // If nothing is present, don't write out anything
 972         if (cSubList.empty())
 973         {
 974             return;
 975         }
 976         else
 977         {
 978             // Take over contents
 979             entityList = cSubList;
 980         }
 981     }
 982     else
 983     {
 984         // Conventional output
 985         entityList = cList;
 986     }
 987
 988     if (useOldPoints)
 989     {
 990         if (useOldConnectivity)
 991         {
 992             // Use points from polyMesh
 993             meshOps::writeVTK
 994             (
 995                 (*this),
 996                 name,
 997                 entityList,
 998                 primitiveType,
 999                 polyMesh::points(),
1000                 polyMesh::edges(),
1001                 polyMesh::faces(),
1002                 polyMesh::cells(),
1003                 polyMesh::faceOwner(),
1004                 scalField,
1005                 lablField
1006             );
1007         }
1008         else
1009         {
1010             meshOps::writeVTK
1011             (
1012                 (*this),
1013                 name,
1014                 entityList,
1015                 primitiveType,
1016                 oldPoints_,
1017                 edges_,
1018                 faces_,
1019                 cells_,
1020                 owner_,
1021                 scalField,
1022                 lablField
1023             );
1024         }
1025     }
1026     else
1027     {
1028         meshOps::writeVTK
1029         (
1030             (*this),
1031             name,
1032             entityList,
1033             primitiveType,
1034             points_,
1035             edges_,
1036             faces_,
1037             cells_,
1038             owner_,
1039             scalField,
1040             lablField
1041         );
1042     }
1043 }
1044
1045
1046 // Return the status report interval
1047 scalar dynamicTopoFvMesh::reportInterval() const
1048 {
1049     // Default to 1 second
1050     scalar interval = 1.0;
1051
1052     const dictionary& meshSubDict = dict_.subDict("dynamicTopoFvMesh");
1053
1054     if (meshSubDict.found("reportInterval") || mandatory_)
1055     {
1056         interval = readScalar(meshSubDict.lookup("reportInterval"));
1057
1058         // Prevent reports if necessary
1059         if (interval < VSMALL)
1060         {
1061             interval = GREAT;
1062         }
1063     }
1064
1065     return interval;
1066 }
1067
1068
1069 // Check the state of connectivity lists
1070 void dynamicTopoFvMesh::checkConnectivity(const label maxErrors) const
1071 {
1072     label nFailedChecks = 0;
1073
1074     messageStream ConnectivityWarning
1075     (
1076         "dynamicTopoFvMesh Connectivity Warning",
1077         messageStream::SERIOUS,
1078         maxErrors
1079     );
1080
1081     // Check face-label ranges
1082     Pout<< "Checking index ranges...";
1083
1084     forAll(edges_, edgeI)
1085     {
1086         const edge& curEdge = edges_[edgeI];
1087
1088         if (curEdge == edge(-1, -1))
1089         {
1090             continue;
1091         }
1092
1093         if
1094         (
1095             curEdge[0] < 0 || curEdge[0] > (points_.size()-1) ||
1096             curEdge[1] < 0 || curEdge[1] > (points_.size()-1)
1097         )
1098         {
1099             Pout<< "Edge " << edgeI
1100                 << " contains vertex labels out of range: "
1101                 << curEdge
1102                 << " Max point index = " << (points_.size()-1) << endl;
1103
1104             nFailedChecks++;
1105
1106             ConnectivityWarning()
1107                 << nl << "Edge-point connectivity is inconsistent."
1108                 << endl;
1109         }
1110
1111         // Check for unique point-labels
1112         if (curEdge[0] == curEdge[1])
1113         {
1114             Pout<< "Edge " << edgeI
1115                 << " contains identical vertex labels: "
1116                 << curEdge
1117                 << endl;
1118
1119             nFailedChecks++;
1120
1121             ConnectivityWarning()
1122                 << nl << "Edge-point connectivity is inconsistent."
1123                 << endl;
1124         }
1125     }
1126
1127     label allPoints = points_.size();
1128     labelList nPointFaces(allPoints, 0);
1129
1130     forAll(faces_, faceI)
1131     {
1132         const face& curFace = faces_[faceI];
1133
1134         if (curFace.empty())
1135         {
1136             continue;
1137         }
1138
1139         if (min(curFace) < 0 || max(curFace) > (points_.size()-1))
1140         {
1141             Pout<< "Face " << faceI
1142                 << " contains vertex labels out of range: "
1143                 << curFace
1144                 << " Max point index = " << (points_.size()-1) << endl;
1145
1146             nFailedChecks++;
1147
1148             ConnectivityWarning()
1149                 << nl << "Face-point connectivity is inconsistent."
1150                 << endl;
1151         }
1152
1153         // Check for unique point-labels
1154         labelHashSet uniquePoints;
1155
1156         forAll(curFace, pointI)
1157         {
1158             bool inserted = uniquePoints.insert(curFace[pointI]);
1159
1160             if (!inserted)
1161             {
1162                 Pout<< "Face " << faceI
1163                     << " contains identical vertex labels: "
1164                     << curFace
1165                     << endl;
1166
1167                 nFailedChecks++;
1168
1169                 ConnectivityWarning()
1170                     << nl << "Face-point connectivity is inconsistent."
1171                     << endl;
1172             }
1173         }
1174
1175         // Count faces per point
1176         forAll(curFace, pointI)
1177         {
1178             nPointFaces[curFace[pointI]]++;
1179         }
1180
1181         // Ensure that cells on either side of this face
1182         // share just one face.
1183         if (neighbour_[faceI] > -1)
1184         {
1185             const cell& ownCell = cells_[owner_[faceI]];
1186             const cell& neiCell = cells_[neighbour_[faceI]];
1187
1188             label nCommon = 0;
1189
1190             forAll(ownCell, fi)
1191             {
1192                 if (findIndex(neiCell, ownCell[fi]) > -1)
1193                 {
1194                     nCommon++;
1195                 }
1196             }
1197
1198             if (nCommon != 1)
1199             {
1200                 Pout<< "Cells: " << nl
1201                     << '\t' << owner_[faceI] << ":: " << ownCell << nl
1202                     << '\t' << neighbour_[faceI] << " :: " << neiCell << nl
1203                     << " share multiple faces. "
1204                     << endl;
1205
1206                 nFailedChecks++;
1207
1208                 ConnectivityWarning()
1209                     << nl << "Cell-Face connectivity is inconsistent."
1210                     << endl;
1211             }
1212         }
1213     }
1214
1215     label allFaces = faces_.size();
1216     labelList nCellsPerFace(allFaces, 0);
1217
1218     forAll(cells_, cellI)
1219     {
1220         const cell& curCell = cells_[cellI];
1221
1222         if (curCell.empty())
1223         {
1224             continue;
1225         }
1226
1227         if (min(curCell) < 0 || max(curCell) > (faces_.size()-1))
1228         {
1229             Pout<< "Cell " << cellI
1230                 << " contains face labels out of range: " << curCell
1231                 << " Max face index = " << (faces_.size()-1) << endl;
1232
1233             nFailedChecks++;
1234
1235             ConnectivityWarning()
1236                 << nl << "Cell-Face connectivity is inconsistent."
1237                 << endl;
1238         }
1239
1240         // Check for unique face-labels
1241         labelHashSet uniqueFaces;
1242
1243         forAll(curCell, faceI)
1244         {
1245             bool inserted = uniqueFaces.insert(curCell[faceI]);
1246
1247             if (!inserted)
1248             {
1249                 Pout<< "Cell " << cellI
1250                     << " contains identical face labels: "
1251                     << curCell
1252                     << endl;
1253
1254                 nFailedChecks++;
1255
1256                 ConnectivityWarning()
1257                     << nl << "Cell-Face connectivity is inconsistent."
1258                     << endl;
1259             }
1260
1261             // Count cells per face
1262             nCellsPerFace[curCell[faceI]]++;
1263         }
1264     }
1265
1266     Pout<< "Done." << endl;
1267
1268     Pout<< "Checking face-cell connectivity...";
1269
1270     forAll(nCellsPerFace, faceI)
1271     {
1272         // This might be a deleted face
1273         if (faceI < nOldFaces_)
1274         {
1275             if (reverseFaceMap_[faceI] == -1)
1276             {
1277                 continue;
1278             }
1279         }
1280         else
1281         if (deletedFaces_.found(faceI))
1282         {
1283             continue;
1284         }
1285
1286         // Determine patch
1287         label uPatch = whichPatch(faceI);
1288
1289         if (nCellsPerFace[faceI] == 0)
1290         {
1291             // Looks like this is really an unused face.
1292             Pout<< "Face " << faceI << " :: " << faces_[faceI]
1293                 << " is unused. Patch: "
1294                 << (uPatch > -1 ? boundaryMesh()[uPatch].name() : "Internal")
1295                 << nl << endl;
1296
1297             nFailedChecks++;
1298
1299             ConnectivityWarning()
1300                 << nl << "Cell-Face connectivity is inconsistent."
1301                 << endl;
1302         }
1303         else
1304         if (nCellsPerFace[faceI] != 2 && uPatch == -1)
1305         {
1306             // Internal face is not shared by exactly two cells
1307             Pout<< "Internal Face " << faceI
1308                 << " :: " << faces_[faceI]
1309                 << " Owner: " << owner_[faceI]
1310                 << " Neighbour: " << neighbour_[faceI]
1311                 << " is multiply connected." << nl
1312                 << " nCellsPerFace: " << nCellsPerFace[faceI] << nl
1313                 << " Patch: Internal" << nl
1314                 << endl;
1315
1316             // Loop through cells and find another instance
1317             forAll(cells_, cellI)
1318             {
1319                 if (findIndex(cells_[cellI], faceI) > -1)
1320                 {
1321                     Pout<< "  Cell: " << cellI
1322                         << "  :: " << cells_[cellI]
1323                         << endl;
1324                 }
1325             }
1326
1327             nFailedChecks++;
1328
1329             ConnectivityWarning()
1330                 << nl << "Cell-Face connectivity is inconsistent."
1331                 << endl;
1332         }
1333         else
1334         if (nCellsPerFace[faceI] != 1 && uPatch > -1)
1335         {
1336             // Boundary face is not shared by exactly one cell
1337             Pout<< "Boundary Face " << faceI
1338                 << " :: " << faces_[faceI]
1339                 << " Owner: " << owner_[faceI]
1340                 << " Neighbour: " << neighbour_[faceI]
1341                 << " is multiply connected." << nl
1342                 << " nCellsPerFace: " << nCellsPerFace[faceI] << nl
1343                 << " Patch: " << boundaryMesh()[uPatch].name() << nl
1344                 << endl;
1345
1346             // Loop through cells and find another instance
1347             forAll(cells_, cellI)
1348             {
1349                 if (findIndex(cells_[cellI], faceI) > -1)
1350                 {
1351                     Pout<< "  Cell: " << cellI
1352                         << "  :: " << cells_[cellI]
1353                         << endl;
1354                 }
1355             }
1356
1357             nFailedChecks++;
1358
1359             ConnectivityWarning()
1360                 << nl << "Cell-Face connectivity is inconsistent."
1361                 << endl;
1362         }
1363     }
1364
1365     Pout<< "Done." << endl;
1366
1367     Pout<< "Checking for unused points...";
1368
1369     forAll(nPointFaces, pointI)
1370     {
1371         if (nPointFaces[pointI] == 0)
1372         {
1373             // This might be a deleted point.
1374             if (pointI < nOldPoints_)
1375             {
1376                 if (reversePointMap_[pointI] == -1)
1377                 {
1378                     continue;
1379                 }
1380             }
1381             else
1382             {
1383                 if (deletedPoints_.found(pointI))
1384                 {
1385                     continue;
1386                 }
1387             }
1388
1389             // Looks like this is really an unused point.
1390             Pout<< nl << nl << "Point " << pointI
1391                 << " is unused. " << endl;
1392
1393             nFailedChecks++;
1394
1395             ConnectivityWarning()
1396                 << nl << "Point-Face connectivity is inconsistent."
1397                 << endl;
1398         }
1399     }
1400
1401     Pout<< "Done." << endl;
1402
1403     Pout<< "Checking edge-face connectivity...";
1404
1405     label allEdges = edges_.size();
1406     labelList nEdgeFaces(allEdges, 0);
1407
1408     forAll(faceEdges_, faceI)
1409     {
1410         const labelList& faceEdges = faceEdges_[faceI];
1411
1412         if (faceEdges.empty())
1413         {
1414             continue;
1415         }
1416
1417         // Check consistency of face-edge-points as well
1418         edgeList eList = faces_[faceI].edges();
1419
1420         forAll(faceEdges,edgeI)
1421         {
1422             nEdgeFaces[faceEdges[edgeI]]++;
1423
1424             // Check if this edge actually belongs to this face
1425             bool found = false;
1426             const edge& edgeToCheck = edges_[faceEdges[edgeI]];
1427
1428             forAll(eList, edgeII)
1429             {
1430                 if (edgeToCheck == eList[edgeII])
1431                 {
1432                     found = true;
1433                     break;
1434                 }
1435             }
1436
1437             if (!found)
1438             {
1439                 Pout<< nl << nl << "Edge: " << faceEdges[edgeI]
1440                     << ": " << edgeToCheck << nl
1441                     << "was not found in face: " << faceI
1442                     << ": " << faces_[faceI] << nl
1443                     << "faceEdges: " << faceEdges
1444                     << endl;
1445
1446                 nFailedChecks++;
1447
1448                 ConnectivityWarning()
1449                     << nl << "Edge-Face connectivity is inconsistent."
1450                     << endl;
1451             }
1452         }
1453     }
1454
1455     label nInternalEdges = 0;
1456     DynamicList<label> bPatchIDs(10);
1457     labelList patchInfo(boundaryMesh().size(), 0);
1458
1459     forAll(edgeFaces_, edgeI)
1460     {
1461         const labelList& edgeFaces = edgeFaces_[edgeI];
1462
1463         if (edgeFaces.empty())
1464         {
1465             continue;
1466         }
1467
1468         if (edgeFaces.size() != nEdgeFaces[edgeI])
1469         {
1470             Pout<< nl << nl << "Edge: " << edgeI << ": " << edges_[edgeI]
1471                 << ": edgeFaces: " << edgeFaces
1472                 << nl << UIndirectList<face>(faces_, edgeFaces)
1473                 << nl << " Expected nFaces: " << nEdgeFaces[edgeI]
1474                 << endl;
1475
1476             nFailedChecks++;
1477
1478             ConnectivityWarning()
1479                 << nl << "Edge-Face connectivity is inconsistent."
1480                 << endl;
1481         }
1482
1483         label nBF = 0;
1484         bPatchIDs.clear();
1485
1486         // Check if this edge belongs to faceEdges for each face
1487         forAll(edgeFaces, faceI)
1488         {
1489             const labelList& faceEdges = faceEdges_[edgeFaces[faceI]];
1490
1491             if (findIndex(faceEdges, edgeI) == -1)
1492             {
1493                 Pout<< nl << nl << "Edge: " << edgeI << ": " << edges_[edgeI]
1494                     << ", edgeFaces: " << edgeFaces
1495                     << nl << UIndirectList<face>(faces_, edgeFaces)
1496                     << "was not found in faceEdges of face: "
1497                     << edgeFaces[faceI] << ": " << faces_[edgeFaces[faceI]]
1498                     << nl << "faceEdges: " << faceEdges
1499                     << nl << UIndirectList<edge>(edges_, faceEdges)
1500                     << endl;
1501
1502                 nFailedChecks++;
1503
1504                 ConnectivityWarning()
1505                     << nl << "Edge-Face connectivity is inconsistent."
1506                     << endl;
1507             }
1508
1509             if (neighbour_[edgeFaces[faceI]] == -1)
1510             {
1511                 // Add to list of patch IDs
1512                 bPatchIDs.append(whichPatch(edgeFaces[faceI]));
1513
1514                 nBF++;
1515             }
1516         }
1517
1518         if (nBF == 0)
1519         {
1520             nInternalEdges++;
1521
1522             // Check if this edge is actually internal.
1523             if (whichEdgePatch(edgeI) >= 0)
1524             {
1525                 Pout<< "Edge: " << edgeI
1526                     << ": " << edges_[edgeI] << " is internal, "
1527                     << " but patch is specified as: "
1528                     << whichEdgePatch(edgeI)
1529                     << endl;
1530
1531                 nFailedChecks++;
1532             }
1533         }
1534         else
1535         {
1536             label patchID = whichEdgePatch(edgeI);
1537
1538             // Check if this edge is actually on a boundary.
1539             if (patchID < 0)
1540             {
1541                 Pout<< "Edge: " << edgeI
1542                     << ": " << edges_[edgeI]
1543                     << " is on a boundary, but patch is specified as: "
1544                     << patchID << endl;
1545
1546                 nFailedChecks++;
1547             }
1548             else
1549             {
1550                 patchInfo[patchID]++;
1551             }
1552
1553             bool failedManifoldCheck = false;
1554
1555             if (nBF > 2)
1556             {
1557                 if (Pstream::parRun())
1558                 {
1559                     // Pinched manifolds should be allowed in parallel
1560                     failedManifoldCheck = false;
1561                 }
1562                 else
1563                 {
1564                     failedManifoldCheck = true;
1565                 }
1566             }
1567
1568             if (failedManifoldCheck)
1569             {
1570                 // Write out for post-processing
1571                 forAll(bPatchIDs, faceI)
1572                 {
1573                     if (bPatchIDs[faceI] == -1)
1574                     {
1575                         Pout<< " Edge: " << edgeI
1576                             << " Face Patch: Internal" << nl;
1577                     }
1578                     else
1579                     {
1580                         Pout<< " Edge: " << edgeI
1581                             << " Face Patch: "
1582                             << boundaryMesh()[bPatchIDs[faceI]].name() << nl;
1583                     }
1584                 }
1585
1586                 Pout<< endl;
1587
1588                 writeVTK("pinched_" + Foam::name(edgeI), edgeFaces, 2);
1589
1590                 Pout<< "Edge: " << edgeI
1591                     << ": " << edges_[edgeI]
1592                     << " has " << nBF
1593                     << " boundary faces connected to it." << nl
1594                     << " Pinched manifolds are not allowed."
1595                     << endl;
1596
1597                 nFailedChecks++;
1598             }
1599         }
1600     }
1601
1602     if (nInternalEdges != nInternalEdges_)
1603     {
1604         Pout<< nl << "Internal edge-count is inconsistent." << nl
1605             << " Counted internal edges: " << nInternalEdges
1606             << " Actual count: " << nInternalEdges_ << endl;
1607
1608         nFailedChecks++;
1609     }
1610
1611     forAll(patchInfo, patchI)
1612     {
1613         if (patchInfo[patchI] != edgePatchSizes_[patchI])
1614         {
1615             Pout<< "Patch-count is inconsistent." << nl
1616                 << " Patch: " << patchI
1617                 << " Counted edges: " << patchInfo[patchI]
1618                 << " Actual count: " << edgePatchSizes_[patchI] << endl;
1619
1620             nFailedChecks++;
1621         }
1622     }
1623
1624     // Check added edge patches to ensure that it is consistent
1625     forAllConstIter(Map<label>, addedEdgePatches_, aepIter)
1626     {
1627         label key = aepIter.key();
1628         label patch = aepIter();
1629
1630         label nBF = 0;
1631         const labelList& edgeFaces = edgeFaces_[key];
1632
1633         // Check if any faces on boundaries
1634         forAll(edgeFaces, faceI)
1635         {
1636             if (neighbour_[edgeFaces[faceI]] == -1)
1637             {
1638                 nBF++;
1639             }
1640         }
1641
1642         if ((patch < 0) && (nBF > 0))
1643         {
1644             Pout<< nl << nl << "Edge: " << key
1645                 << "::" << edges_[key]
1646                 << ", edgeFaces: " << edgeFaces
1647                 << " is internal, but contains boundary faces."
1648                 << endl;
1649
1650             nFailedChecks++;
1651         }
1652
1653         if ((patch >= 0) && (nBF != 2))
1654         {
1655             if (!Pstream::parRun())
1656             {
1657                 Pout<< nl << nl << "Edge: " << key
1658                     << "::" << edges_[key]
1659                     << ", edgeFaces: " << edgeFaces
1660                     << " is on a boundary patch, but doesn't contain"
1661                     << " two boundary faces."
1662                     << endl;
1663
1664                 nFailedChecks++;
1665             }
1666         }
1667     }
1668
1669     Pout<< "Done." << endl;
1670
1671     // Check coupled-patch sizes
1672     forAll(patchCoupling_, patchI)
1673     {
1674         if (patchCoupling_(patchI))
1675         {
1676             const coupleMap& cMap = patchCoupling_[patchI].map();
1677
1678             label mSize = patchSizes_[cMap.masterIndex()];
1679             label sSize = patchSizes_[cMap.slaveIndex()];
1680
1681             if (mSize != sSize)
1682             {
1683                 Pout<< "Coupled patch-count is inconsistent." << nl
1684                     << " Master Patch: " << cMap.masterIndex()
1685                     << " Count: " << mSize << nl
1686                     << " Slave Patch: " << cMap.slaveIndex()
1687                     << " Count: " << sSize
1688                     << endl;
1689
1690                 nFailedChecks++;
1691             }
1692         }
1693     }
1694
1695     if (!twoDMesh_)
1696     {
1697         Pout<< "Checking point-edge connectivity...";
1698
1699         label allPoints = points_.size();
1700         List<labelHashSet> hlPointEdges(allPoints);
1701
1702         forAll(edges_, edgeI)
1703         {
1704             if (edgeFaces_[edgeI].size())
1705             {
1706                 hlPointEdges[edges_[edgeI][0]].insert(edgeI);
1707                 hlPointEdges[edges_[edgeI][1]].insert(edgeI);
1708             }
1709         }
1710
1711         forAll(pointEdges_, pointI)
1712         {
1713             const labelList& pointEdges = pointEdges_[pointI];
1714
1715             if (pointEdges.empty())
1716             {
1717                 continue;
1718             }
1719
1720             forAll(pointEdges, edgeI)
1721             {
1722                 if (!hlPointEdges[pointI].found(pointEdges[edgeI]))
1723                 {
1724                     Pout<< nl << nl << "Point: " << pointI << nl
1725                         << "pointEdges: " << pointEdges << nl
1726                         << "hlPointEdges: " << hlPointEdges[pointI]
1727                         << endl;
1728
1729                     nFailedChecks++;
1730
1731                     ConnectivityWarning()
1732                         << nl << "Point-Edge connectivity is inconsistent."
1733                         << endl;
1734                 }
1735             }
1736
1737             // Do a size check as well
1738             if
1739             (
1740                 hlPointEdges[pointI].size() != pointEdges.size() ||
1741                 pointEdges.size() == 1
1742             )
1743             {
1744                 Pout<< nl << nl << "Point: " << pointI << nl
1745                     << "pointEdges: " << pointEdges << nl
1746                     << "hlPointEdges: " << hlPointEdges[pointI]
1747                     << endl;
1748
1749                 nFailedChecks++;
1750
1751                 ConnectivityWarning()
1752                     << nl << "Size inconsistency."
1753                     << nl << "Point-Edge connectivity is inconsistent."
1754                     << endl;
1755             }
1756         }
1757
1758         Pout<< "Done." << endl;
1759     }
1760
1761     Pout<< "Checking cell-point connectivity...";
1762
1763     // Loop through all cells and construct cell-to-node
1764     label cIndex = 0;
1765     label allCells = cells_.size();
1766     labelList cellIndex(allCells);
1767     List<labelHashSet> cellToNode(allCells);
1768
1769     forAll(cells_, cellI)
1770     {
1771         const cell& thisCell = cells_[cellI];
1772
1773         if (thisCell.empty())
1774         {
1775             continue;
1776         }
1777
1778         cellIndex[cIndex] = cellI;
1779
1780         forAll(thisCell, faceI)
1781         {
1782             const labelList& fEdges = faceEdges_[thisCell[faceI]];
1783
1784             forAll(fEdges, edgeI)
1785             {
1786                 const edge& thisEdge = edges_[fEdges[edgeI]];
1787
1788                 if (!cellToNode[cIndex].found(thisEdge[0]))
1789                 {
1790                     cellToNode[cIndex].insert(thisEdge[0]);
1791                 }
1792
1793                 if (!cellToNode[cIndex].found(thisEdge[1]))
1794                 {
1795                     cellToNode[cIndex].insert(thisEdge[1]);
1796                 }
1797             }
1798         }
1799
1800         cIndex++;
1801     }
1802
1803     // Resize the lists
1804     cellIndex.setSize(cIndex);
1805     cellToNode.setSize(cIndex);
1806
1807     // Preliminary check for size
1808     forAll(cellToNode, cellI)
1809     {
1810         // Check for hexahedral cells
1811         if
1812         (
1813             (cellToNode[cellI].size() == 8) &&
1814             (cells_[cellIndex[cellI]].size() == 6)
1815         )
1816         {
1817             continue;
1818         }
1819
1820         if
1821         (
1822             (cellToNode[cellI].size() != 6 && twoDMesh_) ||
1823             (cellToNode[cellI].size() != 4 && !twoDMesh_)
1824         )
1825         {
1826             Pout<< nl << "Warning: Cell: "
1827                 << cellIndex[cellI] << " is inconsistent. "
1828                 << endl;
1829
1830             const cell& failedCell = cells_[cellIndex[cellI]];
1831
1832             Pout<< "Cell faces: " << failedCell << endl;
1833
1834             forAll(failedCell, faceI)
1835             {
1836                 Pout<< "\tFace: " << failedCell[faceI]
1837                     << " :: " << faces_[failedCell[faceI]]
1838                     << endl;
1839
1840                 const labelList& fEdges = faceEdges_[failedCell[faceI]];
1841
1842                 forAll(fEdges, edgeI)
1843                 {
1844                     Pout<< "\t\tEdge: " << fEdges[edgeI]
1845                         << " :: " << edges_[fEdges[edgeI]]
1846                         << endl;
1847                 }
1848             }
1849
1850             nFailedChecks++;
1851         }
1852     }
1853
1854     Pout<< "Done." << endl;
1855
1856     if (nFailedChecks)
1857     {
1858         FatalErrorIn
1859         (
1860             "void dynamicTopoFvMesh::checkConnectivity"
1861             "(const label maxErrors) const"
1862         )
1863             << nFailedChecks << " failures were found in connectivity."
1864             << abort(FatalError);
1865     }
1866 }
1867
1868
1869 // Utility method to check the quality
1870 // of a triangular face after bisection.
1871 //  - Returns 'true' if the bisection in NOT feasible.
1872 bool dynamicTopoFvMesh::checkBisection
1873 (
1874     const label fIndex,
1875     const label bFaceIndex,
1876     bool forceOp
1877 ) const
1878 {
1879     scalar bisectionQuality = GREAT, minArea = GREAT;
1880
1881     label commonEdge = -1;
1882
1883     // Find the common edge index
1884     meshOps::findCommonEdge
1885     (
1886         bFaceIndex,
1887         fIndex,
1888         faceEdges_,
1889         commonEdge
1890     );
1891
1892     // Fetch the edge
1893     const edge& checkEdge = edges_[commonEdge];
1894     const face& checkFace = faces_[bFaceIndex];
1895
1896     // Compute old / new mid-points
1897     point mpOld =
1898     (
1899         linePointRef
1900         (
1901             oldPoints_[checkEdge.start()],
1902             oldPoints_[checkEdge.end()]
1903         ).centre()
1904     );
1905
1906     point mpNew =
1907     (
1908         linePointRef
1909         (
1910             points_[checkEdge.start()],
1911             points_[checkEdge.end()]
1912         ).centre()
1913     );
1914
1915     // Find the isolated point on the face
1916     label iPoint = -1, nPoint = -1;
1917
1918     meshOps::findIsolatedPoint
1919     (
1920         checkFace,
1921         checkEdge,
1922         iPoint,
1923         nPoint
1924     );
1925
1926     // Find the other point
1927     label oPoint =
1928     (
1929         (nPoint == checkEdge.start()) ?
1930         checkEdge.end() : checkEdge.start()
1931     );
1932
1933     // Configure old / new triangle faces
1934     FixedList<FixedList<point, 3>, 2> tfNew, tfOld;
1935
1936     tfNew[0][0] = mpNew;
1937     tfNew[0][1] = points_[iPoint];
1938     tfNew[0][2] = points_[nPoint];
1939
1940     tfOld[0][0] = mpOld;
1941     tfOld[0][1] = oldPoints_[iPoint];
1942     tfOld[0][2] = oldPoints_[nPoint];
1943
1944     tfNew[1][0] = points_[oPoint];
1945     tfNew[1][1] = points_[iPoint];
1946     tfNew[1][2] = mpNew;
1947
1948     tfOld[1][0] = oldPoints_[oPoint];
1949     tfOld[1][1] = oldPoints_[iPoint];
1950     tfOld[1][2] = mpOld;
1951
1952     // Assume XY plane here
1953     vector n = vector(0,0,1);
1954
1955     forAll(tfNew, fI)
1956     {
1957         // Configure triangles
1958         triPointRef tprNew(tfNew[fI][0], tfNew[fI][1], tfNew[fI][2]);
1959         triPointRef tprOld(tfOld[fI][0], tfOld[fI][1], tfOld[fI][2]);
1960
1961         scalar tQuality =
1962         (
1963             tprNew.quality() *
1964             (
1965                 Foam::sign
1966                 (
1967                     tprNew.normal() &
1968                     ((tprNew.centre() & n) * n)
1969                 )
1970             )
1971         );
1972
1973         scalar oldArea =
1974         (
1975             mag(tprOld.normal()) *
1976             (
1977                 Foam::sign
1978                 (
1979                     tprOld.normal() &
1980                     ((tprOld.centre() & n) * n)
1981                 )
1982             )
1983         );
1984
1985         // Update statistics
1986         minArea = Foam::min(minArea, oldArea);
1987         bisectionQuality = Foam::min(bisectionQuality, tQuality);
1988     }
1989
1990     // Final quality check
1991     if (bisectionQuality < sliverThreshold_ && !forceOp)
1992     {
1993         return true;
1994     }
1995
1996     // Negative quality is a no-no
1997     if (bisectionQuality < 0.0)
1998     {
1999         return true;
2000     }
2001
2002     // Negative old-area is also a no-no
2003     if (minArea < 0.0)
2004     {
2005         return true;
2006     }
2007
2008     // No problems, so a bisection is feasible.
2009     return false;
2010 }
2011
2012
2013 // Utility method to check whether the faces in triFaces will yield
2014 // valid triangles when 'pointIndex' is moved to 'newPoint'.
2015 //  - The routine performs metric-based checks.
2016 //  - Returns 'true' if the collapse in NOT feasible.
2017 //  - Does not reference member data, because this function
2018 //    is also used on subMeshes
2019 bool dynamicTopoFvMesh::checkCollapse
2020 (
2021     const dynamicTopoFvMesh& mesh,
2022     const labelList& triFaces,
2023     const FixedList<label,2>& c0BdyIndex,
2024     const FixedList<label,2>& c1BdyIndex,
2025     const FixedList<label,2>& pointIndex,
2026     const FixedList<point,2>& newPoint,
2027     const FixedList<point,2>& oldPoint,
2028     scalar& collapseQuality,
2029     const bool checkNeighbour,
2030     bool forceOp
2031 )
2032 {
2033     // Reset input
2034     collapseQuality = GREAT;
2035     scalar minArea = GREAT;
2036
2037     forAll(triFaces, indexI)
2038     {
2039         if
2040         (
2041             (triFaces[indexI] == c0BdyIndex[0])
2042          || (triFaces[indexI] == c0BdyIndex[1])
2043         )
2044         {
2045             continue;
2046         }
2047
2048         if (checkNeighbour)
2049         {
2050             if
2051             (
2052                 (triFaces[indexI] == c1BdyIndex[0])
2053              || (triFaces[indexI] == c1BdyIndex[1])
2054             )
2055             {
2056                 continue;
2057             }
2058         }
2059
2060         const face& checkFace = mesh.faces_[triFaces[indexI]];
2061
2062         // Configure a triangle face
2063         FixedList<point, 3> tFNew(vector::zero);
2064         FixedList<point, 3> tFOld(vector::zero);
2065
2066         // Make necessary replacements
2067         forAll(checkFace, pointI)
2068         {
2069             tFNew[pointI] = mesh.points_[checkFace[pointI]];
2070             tFOld[pointI] = mesh.oldPoints_[checkFace[pointI]];
2071
2072             if (checkFace[pointI] == pointIndex[0])
2073             {
2074                 tFNew[pointI] = newPoint[0];
2075                 tFOld[pointI] = oldPoint[0];
2076             }
2077
2078             if (checkFace[pointI] == pointIndex[1])
2079             {
2080                 tFNew[pointI] = newPoint[1];
2081                 tFOld[pointI] = oldPoint[1];
2082             }
2083         }
2084
2085         // Configure triangles
2086         triPointRef tprNew(tFNew[0], tFNew[1], tFNew[2]);
2087         triPointRef tprOld(tFOld[0], tFOld[1], tFOld[2]);
2088
2089         // Assume XY plane here
2090         vector n = vector(0,0,1);
2091
2092         // Compute the quality.
2093         // Assume centre-plane passes through origin
2094         scalar tQuality =
2095         (
2096             tprNew.quality() *
2097             (
2098                 Foam::sign
2099                 (
2100                     tprNew.normal() &
2101                     ((tprNew.centre() & n) * n)
2102                 )
2103             )
2104         );
2105
2106         scalar oldArea =
2107         (
2108             mag(tprOld.normal()) *
2109             (
2110                 Foam::sign
2111                 (
2112                     tprOld.normal() &
2113                     ((tprOld.centre() & n) * n)
2114                 )
2115             )
2116         );
2117
2118         // Update statistics
2119         minArea = Foam::min(minArea, oldArea);
2120         collapseQuality = Foam::min(collapseQuality, tQuality);
2121     }
2122
2123     // Final quality check
2124     if (collapseQuality < mesh.sliverThreshold_ && !forceOp)
2125     {
2126         if (debug > 3)
2127         {
2128             Pout<< " * * * 2D checkCollapse * * * " << nl
2129                 << " collapseQuality: " << collapseQuality
2130                 << " below threshold: " << mesh.sliverThreshold_
2131                 << endl;
2132         }
2133
2134         return true;
2135     }
2136
2137     // Negative quality is a no-no
2138     if (collapseQuality < 0.0)
2139     {
2140         if (forceOp)
2141         {
2142             Pout<< " * * * 2D checkCollapse * * * " << nl
2143                 << " Negative collapseQuality: " << collapseQuality << nl
2144                 << " Operation cannot be forced."
2145                 << endl;
2146         }
2147
2148         return true;
2149     }
2150
2151     // Negative old-area is also a no-no
2152     if (minArea < 0.0)
2153     {
2154         if (forceOp)
2155         {
2156             Pout<< " * * * 2D checkCollapse * * * " << nl
2157                 << " minArea: " << minArea << nl
2158                 << " Operation cannot be forced."
2159                 << endl;
2160         }
2161
2162         return true;
2163     }
2164
2165     // No problems, so a collapse is feasible.
2166     return false;
2167 }
2168
2169
2170 // Utility method to check whether the cell given by 'cellIndex' will yield
2171 // a valid cell when 'pointIndex' is moved to 'newPoint'.
2172 //  - The routine performs metric-based checks.
2173 //  - Returns 'true' if the collapse in NOT feasible, and
2174 //    makes entries in cellsChecked to avoid repetitive checks.
2175 bool dynamicTopoFvMesh::checkCollapse
2176 (
2177     const point& newPoint,
2178     const point& oldPoint,
2179     const label pointIndex,
2180     const label cellIndex,
2181     DynamicList<label>& cellsChecked,
2182     scalar& collapseQuality,
2183     bool forceOp
2184 ) const
2185 {
2186     label faceIndex = -1;
2187     scalar cQuality = 0;
2188     scalar oldVolume = 0;
2189     scalar newVolume = 0;
2190     const cell& cellToCheck = cells_[cellIndex];
2191
2192     // Look for a face that doesn't contain 'pointIndex'
2193     forAll(cellToCheck, faceI)
2194     {
2195         const face& currFace = faces_[cellToCheck[faceI]];
2196
2197         if (currFace.which(pointIndex) < 0)
2198         {
2199             faceIndex = cellToCheck[faceI];
2200             break;
2201         }
2202     }
2203
2204     // Compute cell-volume
2205     const face& faceToCheck = faces_[faceIndex];
2206
2207     if (owner_[faceIndex] == cellIndex)
2208     {
2209         cQuality =
2210         (
2211             tetMetric_
2212             (
2213                 points_[faceToCheck[2]],
2214                 points_[faceToCheck[1]],
2215                 points_[faceToCheck[0]],
2216                 newPoint
2217             )
2218         );
2219
2220         newVolume =
2221         (
2222             tetPointRef
2223             (
2224                 points_[faceToCheck[2]],
2225                 points_[faceToCheck[1]],
2226                 points_[faceToCheck[0]],
2227                 newPoint
2228             ).mag()
2229         );
2230
2231         oldVolume =
2232         (
2233             tetPointRef
2234             (
2235                 oldPoints_[faceToCheck[2]],
2236                 oldPoints_[faceToCheck[1]],
2237                 oldPoints_[faceToCheck[0]],
2238                 oldPoint
2239             ).mag()
2240         );
2241     }
2242     else
2243     {
2244         cQuality =
2245         (
2246             tetMetric_
2247             (
2248                 points_[faceToCheck[0]],
2249                 points_[faceToCheck[1]],
2250                 points_[faceToCheck[2]],
2251                 newPoint
2252             )
2253         );
2254
2255         newVolume =
2256         (
2257             tetPointRef
2258             (
2259                 points_[faceToCheck[0]],
2260                 points_[faceToCheck[1]],
2261                 points_[faceToCheck[2]],
2262                 newPoint
2263             ).mag()
2264         );
2265
2266         oldVolume =
2267         (
2268             tetPointRef
2269             (
2270                 oldPoints_[faceToCheck[0]],
2271                 oldPoints_[faceToCheck[1]],
2272                 oldPoints_[faceToCheck[2]],
2273                 oldPoint
2274             ).mag()
2275         );
2276     }
2277
2278     // Final quality check
2279     if (cQuality < sliverThreshold_ && !forceOp)
2280     {
2281         if (debug > 4)
2282         {
2283             Pout<< " * * * 3D checkCollapse * * * " << nl
2284                 << "\nCollapsing cell: " << cellIndex
2285                 << " containing points:\n"
2286                 << faceToCheck[0] << "," << faceToCheck[1] << ","
2287                 << faceToCheck[2] << "," << pointIndex << nl
2288                 << "will yield a quality of: " << cQuality
2289                 << ", when " << pointIndex
2290                 << " is moved to location: " << nl
2291                 << newPoint
2292                 << endl;
2293         }
2294
2295         return true;
2296     }
2297
2298     // Negative quality is a no-no
2299     if (cQuality < 0.0)
2300     {
2301         if (forceOp)
2302         {
2303             Pout<< " * * * 3D checkCollapse * * * " << nl
2304                 << "\nCollapsing cell: " << cellIndex
2305                 << " containing points:\n"
2306                 << faceToCheck[0] << "," << faceToCheck[1] << ","
2307                 << faceToCheck[2] << "," << pointIndex << nl
2308                 << "will yield a quality of: " << cQuality
2309                 << ", when " << pointIndex
2310                 << " is moved to location: " << nl
2311                 << newPoint << nl
2312                 << "Operation cannot be forced."
2313                 << endl;
2314         }
2315
2316         return true;
2317     }
2318
2319     // Negative old-volume is also a no-no
2320     if (oldVolume < 0.0 || (mag(oldVolume) < mag(0.1 * newVolume)))
2321     {
2322         if (forceOp)
2323         {
2324             Pout<< " * * * 3D checkCollapse * * * " << nl
2325                 << "\nCollapsing cell: " << cellIndex
2326                 << " containing points:\n"
2327                 << faceToCheck[0] << "," << faceToCheck[1] << ","
2328                 << faceToCheck[2] << "," << pointIndex << nl
2329                 << "will yield an old-volume of: " << oldVolume
2330                 << ", when " << pointIndex
2331                 << " is moved to location: " << nl
2332                 << oldPoint << nl
2333                 << "newVolume: " << newVolume << nl
2334                 << "Operation cannot be forced."
2335                 << endl;
2336         }
2337
2338         return true;
2339     }
2340
2341     // No problems, so a collapse is feasible
2342     cellsChecked.append(cellIndex);
2343
2344     // Update input quality
2345     collapseQuality = Foam::min(collapseQuality, cQuality);
2346
2347     return false;
2348 }
2349
2350
2351 } // End namespace Foam
2352
2353 // ************************************************************************* //