src/dynamicMesh/dynamicTopoFvMesh/dynamicTopoFvMeshCheck.C

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