src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/meshOps.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     meshOps
  27
  28 Description
  29     Various utility functions that perform mesh-related operations.
  30
  31 Author
  32     Sandeep Menon
  33     University of Massachusetts Amherst
  34     All rights reserved
  35
  36 \*---------------------------------------------------------------------------*/
  37
  38 #include "objectRegistry.H"
  39 #include "Time.H"
  40 #include "meshOps.H"
  41 #include "ListOps.H"
  42 #include "Pstream.H"
  43 #include "triFace.H"
  44 #include "IOmanip.H"
  45 #include "HashSet.H"
  46 #include "polyMesh.H"
  47 #include "triPointRef.H"
  48 #include "tetPointRef.H"
  49
  50 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  51
  52 namespace Foam
  53 {
  54
  55 namespace meshOps
  56 {
  57
  58 // Utility method to build a hull of cells
  59 // connected to the edge (for 2D simplical meshes)
  60 void constructPrismHull
  61 (
  62     const label eIndex,
  63     const UList<face>& faces,
  64     const UList<cell>& cells,
  65     const UList<label>& owner,
  66     const UList<label>& neighbour,
  67     const UList<labelList>& edgeFaces,
  68     labelList& hullTriFaces,
  69     labelList& hullCells
  70 )
  71 {
  72     labelHashSet cellSet, triFaceSet;
  73
  74     // Obtain references
  75     const labelList& eFaces = edgeFaces[eIndex];
  76
  77     // Loop through edgeFaces and add cells
  78     forAll(eFaces, faceI)
  79     {
  80         label c0 = owner[eFaces[faceI]];
  81         label c1 = neighbour[eFaces[faceI]];
  82
  83         if (!cellSet.found(c0))
  84         {
  85             // Add this cell
  86             cellSet.insert(c0);
  87
  88             // Find associated triFaces and add them too
  89             const cell& cC = cells[c0];
  90
  91             forAll(cC, faceJ)
  92             {
  93                 const face& cF = faces[cC[faceJ]];
  94
  95                 if ((cF.size() == 3) && !(triFaceSet.found(cC[faceJ])))
  96                 {
  97                     triFaceSet.insert(cC[faceJ]);
  98                 }
  99             }
 100         }
 101
 102         if (!cellSet.found(c1) && (c1 != -1))
 103         {
 104             // Add this cell
 105             cellSet.insert(c1);
 106
 107             // Find associated triFaces and add them too
 108             const cell& cC = cells[c1];
 109
 110             forAll(cC, faceJ)
 111             {
 112                 const face& cF = faces[cC[faceJ]];
 113
 114                 if ((cF.size() == 3) && !(triFaceSet.found(cC[faceJ])))
 115                 {
 116                     triFaceSet.insert(cC[faceJ]);
 117                 }
 118             }
 119         }
 120     }
 121
 122     // Obtain lists from hashSets
 123     hullCells = cellSet.toc();
 124     hullTriFaces = triFaceSet.toc();
 125 }
 126
 127
 128 // Utility method to build a hull of cells (and faces)
 129 // around an edge (for 3D simplical meshes)
 130 void constructHull
 131 (
 132     const label eIndex,
 133     const UList<face>& faces,
 134     const UList<edge>& edges,
 135     const UList<cell>& cells,
 136     const UList<label>& owner,
 137     const UList<label>& neighbour,
 138     const UList<labelList>& faceEdges,
 139     const UList<labelList>& edgeFaces,
 140     const labelList& vertexHull,
 141     labelList& hullEdges,
 142     labelList& hullFaces,
 143     labelList& hullCells,
 144     labelListList& ringEntities
 145 )
 146 {
 147     // [1] hullEdges is an ordered list of edge-labels around eIndex,
 148     //     but not connected to it.
 149     //      - Ordering is in the same manner as vertexHull.
 150     // [2] hullFaces is an ordered list of face-labels connected to eIndex.
 151     //      - Ordering is in the same manner as vertexHull.
 152     // [3] hullCells is an ordered list of cell-labels connected to eIndex.
 153     //      - For boundary hulls, the last cell label is -1
 154     // [4] ringEntities are edges and faces connected to eIndex[0] and eIndex[1]
 155     //      - ringEntities[0]: edges connected to eIndex[0]
 156     //      - ringEntities[1]: faces connected to eIndex[0]
 157     //      - ringEntities[2]: edges connected to eIndex[1]
 158     //      - ringEntities[3]: faces connected to eIndex[1]
 159
 160     bool found;
 161     label otherPoint = -1, nextPoint = -1;
 162
 163     // Obtain a reference to this edge, and its edgeFaces
 164     const edge& edgeToCheck = edges[eIndex];
 165     const labelList& eFaces = edgeFaces[eIndex];
 166
 167     // Loop through all faces of this edge and add them to hullFaces
 168     forAll(eFaces, faceI)
 169     {
 170         const face& faceToCheck = faces[eFaces[faceI]];
 171
 172         // Find the isolated point on this face,
 173         // and compare it with vertexHull
 174         meshOps::findIsolatedPoint
 175         (
 176             faceToCheck,
 177             edgeToCheck,
 178             otherPoint,
 179             nextPoint
 180         );
 181
 182         found = false;
 183
 184         forAll(vertexHull, indexI)
 185         {
 186             if (vertexHull[indexI] == otherPoint)
 187             {
 188                 // Fill in the position of this face on the hull
 189                 hullFaces[indexI] = eFaces[faceI];
 190
 191                 // Obtain edges connected to top and bottom
 192                 // vertices of edgeToCheck
 193                 const labelList& fEdges = faceEdges[hullFaces[indexI]];
 194
 195                 forAll(fEdges, edgeI)
 196                 {
 197                     if
 198                     (
 199                         edges[fEdges[edgeI]]
 200                      == edge(edgeToCheck[0], otherPoint)
 201                     )
 202                     {
 203                         ringEntities[0][indexI] = fEdges[edgeI];
 204                     }
 205
 206                     if
 207                     (
 208                         edges[fEdges[edgeI]]
 209                      == edge(edgeToCheck[1], otherPoint)
 210                     )
 211                     {
 212                         ringEntities[2][indexI] = fEdges[edgeI];
 213                     }
 214                 }
 215
 216                 // Depending on the orientation of this face,
 217                 // fill in hull cell indices as well
 218                 if (nextPoint == edgeToCheck[0])
 219                 {
 220                     hullCells[indexI] = owner[eFaces[faceI]];
 221                 }
 222                 else
 223                 if (nextPoint == edgeToCheck[1])
 224                 {
 225                     hullCells[indexI] = neighbour[eFaces[faceI]];
 226                 }
 227                 else
 228                 {
 229                     // Something's terribly wrong
 230                     FatalErrorIn("void meshOps::constructHull()")
 231                         << nl << " Failed to construct hull. "
 232                         << nl << " Possibly not a tetrahedral mesh. "
 233                         << abort(FatalError);
 234                 }
 235
 236                 if (hullCells[indexI] != -1)
 237                 {
 238                     label nextI = vertexHull.fcIndex(indexI);
 239                     label nextHullPoint = vertexHull[nextI];
 240                     const cell& currCell = cells[hullCells[indexI]];
 241
 242                     // Look for the ring-faces
 243                     forAll(currCell, faceI)
 244                     {
 245                         const face& cFace = faces[currCell[faceI]];
 246
 247                         // Check if this face contains edgeToCheck[0]
 248                         if
 249                         (
 250                             triFace::compare
 251                             (
 252                                 triFace(cFace),
 253                                 triFace
 254                                 (
 255                                     edgeToCheck[0],
 256                                     otherPoint,
 257                                     nextHullPoint
 258                                 )
 259                             )
 260                         )
 261                         {
 262                             ringEntities[1][indexI] = currCell[faceI];
 263                         }
 264
 265                         // Check if this face contains edgeToCheck[1]
 266                         if
 267                         (
 268                             triFace::compare
 269                             (
 270                                 triFace(cFace),
 271                                 triFace
 272                                 (
 273                                     edgeToCheck[1],
 274                                     nextHullPoint,
 275                                     otherPoint
 276                                 )
 277                             )
 278                         )
 279                         {
 280                             ringEntities[3][indexI] = currCell[faceI];
 281                         }
 282                     }
 283
 284                     // Scan one the faces for the ring-edge
 285                     const labelList& rFaceEdges =
 286                     (
 287                         faceEdges[ringEntities[1][indexI]]
 288                     );
 289
 290                     forAll(rFaceEdges, edgeI)
 291                     {
 292                         if
 293                         (
 294                             edges[rFaceEdges[edgeI]]
 295                          == edge(otherPoint,nextHullPoint)
 296                         )
 297                         {
 298                             hullEdges[indexI] = rFaceEdges[edgeI];
 299                             break;
 300                         }
 301                     }
 302                 }
 303
 304                 // Done with this index. Break out.
 305                 found = true;
 306                 break;
 307             }
 308         }
 309
 310         // Throw an error if the point wasn't found
 311         if (!found)
 312         {
 313             // Something's terribly wrong
 314             FatalErrorIn("void meshOps::constructHull()")
 315                 << " Failed to construct hull. " << nl
 316                 << " edgeFaces connectivity is inconsistent. " << nl
 317                 << " Edge: " << eIndex << ":: " << edgeToCheck << nl
 318                 << " edgeFaces: " << eFaces << nl
 319                 << " vertexHull: " << vertexHull
 320                 << abort(FatalError);
 321         }
 322     }
 323 }
 324
 325
 326 // Given a set of points and edges, find the shortest path
 327 // between the start and end point, using Dijkstra's algorithm.
 328 //  - Takes a Map of points and edges that use those points.
 329 //  - Edge weights are currently edge-lengths, but can easily be adapted.
 330 //  - Returns true if the endPoint was found by the algorithm.
 331 //  - The Map 'pi' returns a preceding point for every point in 'points'.
 332 //
 333 //  Algorithm is inspired by:
 334 //    Renaud Waldura
 335 //    Dijkstra's Shortest Path Algorithm in Java
 336 //    http://renaud.waldura.com/
 337 bool Dijkstra
 338 (
 339     const Map<point>& points,
 340     const Map<edge>& edges,
 341     const label startPoint,
 342     const label endPoint,
 343     Map<label>& pi
 344 )
 345 {
 346     bool foundEndPoint = false;
 347
 348     // Set of unvisited (Q) / visited (S) points and distances (d)
 349     labelHashSet Q, S;
 350     Map<scalar> d;
 351
 352     // Initialize distances to large values
 353     forAllConstIter(Map<point>, points, pIter)
 354     {
 355         d.insert(pIter.key(), GREAT);
 356     }
 357
 358     // Invert edges to make a local pointEdges list
 359     Map<labelList> localPointEdges;
 360
 361     forAllConstIter(Map<edge>, edges, eIter)
 362     {
 363         const edge& edgeToCheck = eIter();
 364
 365         forAll(edgeToCheck, pointI)
 366         {
 367             if (!localPointEdges.found(edgeToCheck[pointI]))
 368             {
 369                 localPointEdges.insert(edgeToCheck[pointI], labelList(0));
 370             }
 371
 372             meshOps::sizeUpList
 373             (
 374                 eIter.key(),
 375                 localPointEdges[edgeToCheck[pointI]]
 376             );
 377         }
 378     }
 379
 380     // Mark the startPoint as having the smallest distance
 381     d[startPoint] = 0.0;
 382
 383     // Add the startPoint to the list of unvisited points
 384     Q.insert(startPoint);
 385
 386     while (Q.size())
 387     {
 388         // Step 1: Find the node with the smallest distance from the start.
 389         labelHashSet::iterator smallest = Q.begin();
 390
 391         for
 392         (
 393             labelHashSet::iterator iter = ++Q.begin();
 394             iter != Q.end();
 395             iter++
 396         )
 397         {
 398             if (d[iter.key()] < d[smallest.key()])
 399             {
 400                 smallest = iter;
 401             }
 402         }
 403
 404         label pointIndex = smallest.key();
 405         scalar smallestDistance = d[pointIndex];
 406
 407         // Move to the visited points list
 408         S.insert(pointIndex);
 409         Q.erase(pointIndex);
 410
 411         // Step 2: Build a list of points adjacent to pointIndex
 412         //         but not in the visited list
 413         DynamicList<label> adjacentPoints(10);
 414
 415         const labelList& pEdges = localPointEdges[pointIndex];
 416
 417         forAll(pEdges, edgeI)
 418         {
 419             const edge& edgeToCheck = edges[pEdges[edgeI]];
 420
 421             label otherPoint = edgeToCheck.otherVertex(pointIndex);
 422
 423             if (!S.found(otherPoint))
 424             {
 425                 adjacentPoints.append(otherPoint);
 426             }
 427         }
 428
 429         // Step 3: Perform distance-based checks for adjacent points
 430         forAll(adjacentPoints, pointI)
 431         {
 432             label adjPoint = adjacentPoints[pointI];
 433
 434             scalar distance =
 435             (
 436                 mag(points[adjPoint] - points[pointIndex])
 437               + smallestDistance
 438             );
 439
 440             // Check if the end-point has been touched.
 441             if (adjPoint == endPoint)
 442             {
 443                 foundEndPoint = true;
 444             }
 445
 446             if (distance < d[adjPoint])
 447             {
 448                 // Update to the shorter distance
 449                 d[adjPoint] = distance;
 450
 451                 // Update the predecessor
 452                 if (pi.found(adjPoint))
 453                 {
 454                     pi[adjPoint] = pointIndex;
 455                 }
 456                 else
 457                 {
 458                     pi.insert(adjPoint, pointIndex);
 459                 }
 460
 461                 // Add to the list of unvisited points
 462                 Q.insert(adjPoint);
 463             }
 464         }
 465     }
 466
 467     return foundEndPoint;
 468 }
 469
 470
 471 // Select a list of elements from connectivity,
 472 // and output to a VTK format
 473 void writeVTK
 474 (
 475     const polyMesh& mesh,
 476     const word& name,
 477     const labelList& cList,
 478     const label primitiveType,
 479     const UList<point>& meshPoints,
 480     const UList<edge>& edges,
 481     const UList<face>& faces,
 482     const UList<cell>& cells,
 483     const UList<label>& owner,
 484     const UList<scalar>& scalField,
 485     const UList<label>& lablField,
 486     const UList<vector>& vectField
 487 )
 488 {
 489     label nTotalCells = 0;
 490     label nPoints = 0, nCells = 0;
 491
 492     // Estimate a size for points and cellPoints
 493     pointField points(6*cList.size());
 494
 495     // Connectivity lists
 496     labelListList cpList(cList.size());
 497
 498     // Create a map for local points
 499     Map<label> pointMap, reversePointMap, reverseCellMap;
 500
 501     forAll(cList, cellI)
 502     {
 503         if (cList[cellI] < 0)
 504         {
 505             continue;
 506         }
 507
 508         bool isPolyhedron = false;
 509
 510         switch (primitiveType)
 511         {
 512             // Are we looking at points?
 513             case 0:
 514             {
 515                 // Size the list
 516                 cpList[nCells].setSize(1);
 517
 518                 cpList[nCells] = cList[cellI];
 519
 520                 nTotalCells++;
 521
 522                 break;
 523             }
 524
 525             // Are we looking at edges?
 526             case 1:
 527             {
 528                 // Size the list
 529                 cpList[nCells].setSize(2);
 530
 531                 const edge& tEdge = edges[cList[cellI]];
 532
 533                 cpList[nCells][0] = tEdge[0];
 534                 cpList[nCells][1] = tEdge[1];
 535
 536                 nTotalCells += 2;
 537
 538                 break;
 539             }
 540
 541             // Are we looking at faces?
 542             case 2:
 543             {
 544                 const face& tFace = faces[cList[cellI]];
 545
 546                 // Size the list and assign
 547                 cpList[nCells] = tFace;
 548
 549                 nTotalCells += tFace.size();
 550
 551                 break;
 552             }
 553
 554             // Are we looking at cells?
 555             case 3:
 556             {
 557                 const cell& tCell = cells[cList[cellI]];
 558
 559                 // Is face information provided?
 560                 if (faces.size())
 561                 {
 562                     if (tCell.size() == 4)
 563                     {
 564                         // Point-ordering for tetrahedra
 565                         const face& baseFace = faces[tCell[0]];
 566                         const face& checkFace = faces[tCell[1]];
 567
 568                         // Size the list
 569                         cpList[nCells].setSize(4);
 570
 571                         // Get the fourth point
 572                         label apexPoint =
 573                         (
 574                             meshOps::findIsolatedPoint(baseFace, checkFace)
 575                         );
 576
 577                         // Something's wrong with connectivity.
 578                         if (apexPoint == -1)
 579                         {
 580                             FatalErrorIn
 581                             (
 582                                 "void meshOps::writeVTK\n"
 583                                 "(\n"
 584                                 "    const polyMesh& mesh,\n"
 585                                 "    const word& name,\n"
 586                                 "    const labelList& cList,\n"
 587                                 "    const label primitiveType,\n"
 588                                 "    const UList<point>& meshPoints,\n"
 589                                 "    const UList<edge>& edges,\n"
 590                                 "    const UList<face>& faces,\n"
 591                                 "    const UList<cell>& cells,\n"
 592                                 "    const UList<label>& owner\n"
 593                                 "    const UList<scalar>& scalField,\n"
 594                                 "    const UList<label>& lablField,\n"
 595                                 "    const UList<vector>& vectField\n"
 596                                 ")\n"
 597                             )
 598                                 << "Cell: " << cList[cellI]
 599                                 << ":: " << tCell
 600                                 << " has inconsistent connectivity."
 601                                 << abort(FatalError);
 602                         }
 603
 604                         // Write-out in order
 605                         label ownCell = owner[tCell[0]];
 606
 607                         if (ownCell == cList[cellI])
 608                         {
 609                             cpList[nCells][0] = baseFace[2];
 610                             cpList[nCells][1] = baseFace[1];
 611                             cpList[nCells][2] = baseFace[0];
 612                             cpList[nCells][3] = apexPoint;
 613                         }
 614                         else
 615                         {
 616                             cpList[nCells][0] = baseFace[0];
 617                             cpList[nCells][1] = baseFace[1];
 618                             cpList[nCells][2] = baseFace[2];
 619                             cpList[nCells][3] = apexPoint;
 620                         }
 621
 622                         nTotalCells += 4;
 623                     }
 624                     else
 625                     if (tCell.size() > 4)
 626                     {
 627                         // Point ordering for polyhedra
 628                         isPolyhedron = true;
 629
 630                         // First obtain the face count
 631                         label npF = 0;
 632
 633                         forAll(tCell, faceI)
 634                         {
 635                             npF += faces[tCell[faceI]].size();
 636                         }
 637
 638                         // Size the list
 639                         cpList[nCells].setSize(tCell.size() + npF + 1);
 640
 641                         // Fill in facePoints
 642                         label nP = 0;
 643
 644                         // Fill in the number of faces
 645                         cpList[nCells][nP++] = tCell.size();
 646
 647                         forAll(tCell, faceI)
 648                         {
 649                             const face& checkFace = faces[tCell[faceI]];
 650
 651                             if (checkFace.empty())
 652                             {
 653                                 FatalErrorIn("void meshOps::writeVTK()")
 654                                     << " Empty face: " << tCell[faceI]
 655                                     << abort(FatalError);
 656                             }
 657
 658                             // First fill in face size
 659                             cpList[nCells][nP++] = checkFace.size();
 660
 661                             // Next fill in points in order
 662                             if (owner[tCell[faceI]] == cList[cellI])
 663                             {
 664                                 forAll(checkFace, pI)
 665                                 {
 666                                     cpList[nCells][nP++] = checkFace[pI];
 667                                 }
 668                             }
 669                             else
 670                             {
 671                                 forAllReverse(checkFace, pI)
 672                                 {
 673                                     cpList[nCells][nP++] = checkFace[pI];
 674                                 }
 675                             }
 676                         }
 677
 678                         nTotalCells += (tCell.size() + npF + 1);
 679                     }
 680                     else
 681                     {
 682                         FatalErrorIn("void meshOps::writeVTK()")
 683                             << " Wrong cell size: " << tCell << nl
 684                             << " Index: " << cList[cellI] << nl
 685                             << abort(FatalError);
 686                     }
 687                 }
 688                 else
 689                 {
 690                     // No face information.
 691                     // Assume cell-to-node information.
 692                     if (tCell.size() == 4)
 693                     {
 694                         // Build a face out of first three points.
 695                         triPointRef tpr
 696                         (
 697                             meshPoints[tCell[0]],
 698                             meshPoints[tCell[1]],
 699                             meshPoints[tCell[2]]
 700                         );
 701
 702                         // Fetch fourth point
 703                         const point& d = meshPoints[tCell[3]];
 704
 705                         scalar dDotN = ((d - tpr.a()) & tpr.normal());
 706
 707                         // Size the list
 708                         cpList[nCells].setSize(4);
 709
 710                         if (dDotN > 0.0)
 711                         {
 712                             // Correct orientation
 713                             cpList[nCells][0] = tCell[0];
 714                             cpList[nCells][1] = tCell[1];
 715                             cpList[nCells][2] = tCell[2];
 716                             cpList[nCells][3] = tCell[3];
 717                         }
 718                         else
 719                         {
 720                             // Flip triangle
 721                             cpList[nCells][0] = tCell[2];
 722                             cpList[nCells][1] = tCell[1];
 723                             cpList[nCells][2] = tCell[0];
 724                             cpList[nCells][3] = tCell[3];
 725                         }
 726
 727                         nTotalCells += 4;
 728                     }
 729                     else
 730                     {
 731                         FatalErrorIn("void meshOps::writeVTK()")
 732                             << " Wrong cell size: " << tCell << nl
 733                             << " Index: " << cList[cellI] << nl
 734                             << abort(FatalError);
 735                     }
 736                 }
 737
 738                 break;
 739             }
 740
 741             default:
 742             {
 743                 FatalErrorIn("void meshOps::writeVTK() const")
 744                     << " Incorrect primitiveType: "
 745                     << primitiveType
 746                     << abort(FatalError);
 747             }
 748         }
 749
 750         // Renumber to local ordering
 751         if (isPolyhedron)
 752         {
 753             // Polyhedra also contain face sizes, so use a different scheme
 754             label pCtr = 0;
 755
 756             // Fetch number of faces
 757             label npF = cpList[nCells][pCtr++];
 758
 759             for (label i = 0; i < npF; i++)
 760             {
 761                 // Fetch number of points
 762                 label npP = cpList[nCells][pCtr++];
 763
 764                 // Read points in order
 765                 for (label j = 0; j < npP; j++)
 766                 {
 767                     // Check if this point was added to the map
 768                     if (!pointMap.found(cpList[nCells][pCtr]))
 769                     {
 770                         // Resize pointField if necessary
 771                         if (nPoints >= points.size())
 772                         {
 773                             points.setSize(2 * nPoints);
 774                         }
 775
 776                         // Point was not found, so add it
 777                         points[nPoints] = meshPoints[cpList[nCells][pCtr]];
 778
 779                         // Update the map
 780                         pointMap.insert(cpList[nCells][pCtr], nPoints);
 781                         reversePointMap.insert(nPoints, cpList[nCells][pCtr]);
 782
 783                         // Increment the number of points
 784                         nPoints++;
 785                     }
 786
 787                     // Renumber it.
 788                     cpList[nCells][pCtr] = pointMap[cpList[nCells][pCtr]];
 789
 790                     // Increment point counter
 791                     pCtr++;
 792                 }
 793             }
 794         }
 795         else
 796         {
 797             forAll(cpList[nCells], pointI)
 798             {
 799                 // Check if this point was added to the map
 800                 if (!pointMap.found(cpList[nCells][pointI]))
 801                 {
 802                     // Resize pointField if necessary
 803                     if (nPoints >= points.size())
 804                     {
 805                         points.setSize(2 * nPoints);
 806                     }
 807
 808                     // Point was not found, so add it
 809                     points[nPoints] = meshPoints[cpList[nCells][pointI]];
 810
 811                     // Update the map
 812                     pointMap.insert(cpList[nCells][pointI], nPoints);
 813                     reversePointMap.insert(nPoints, cpList[nCells][pointI]);
 814
 815                     // Increment the number of points
 816                     nPoints++;
 817                 }
 818
 819                 // Renumber it.
 820                 cpList[nCells][pointI] = pointMap[cpList[nCells][pointI]];
 821             }
 822         }
 823
 824         // Update the cell map.
 825         reverseCellMap.insert(nCells, cList[cellI]);
 826
 827         nCells++;
 828     }
 829
 830     // Finally write it out
 831     meshOps::writeVTK
 832     (
 833         mesh,
 834         name,
 835         nPoints,
 836         nCells,
 837         nTotalCells,
 838         points,
 839         cpList,
 840         primitiveType,
 841         reversePointMap,
 842         reverseCellMap,
 843         scalField,
 844         lablField,
 845         vectField
 846     );
 847 }
 848
 849
 850 // Actual routine to write out the VTK file
 851 void writeVTK
 852 (
 853     const polyMesh& mesh,
 854     const word& name,
 855     const label nPoints,
 856     const label nCells,
 857     const label nTotalCells,
 858     const vectorField& points,
 859     const labelListList& cpList,
 860     const label primitiveType,
 861     const Map<label>& reversePointMap,
 862     const Map<label>& reverseCellMap,
 863     const UList<scalar>& scalField,
 864     const UList<label>& lablField,
 865     const UList<vector>& vectField
 866 )
 867 {
 868     // Make the directory
 869     fileName dirName(mesh.time().path()/"VTK"/mesh.time().timeName());
 870
 871     mkDir(dirName);
 872
 873     // Open stream for output
 874     OFstream file(dirName/name+".vtk");
 875
 876     // Write out the header
 877     file << "# vtk DataFile Version 2.0" << nl
 878          << name << ".vtk" << nl
 879          << "ASCII" << nl
 880          << "DATASET UNSTRUCTURED_GRID" << nl
 881          << "POINTS " << nPoints << " double" << nl;
 882
 883     for (label i = 0; i < nPoints; i++)
 884     {
 885         file << setprecision(15)
 886              << points[i].x() << ' '
 887              << points[i].y() << ' '
 888              << points[i].z() << ' '
 889              << nl;
 890     }
 891
 892     file << "CELLS " << nCells << " " << nTotalCells + nCells << endl;
 893
 894     if (cpList.size())
 895     {
 896         forAll(cpList, i)
 897         {
 898             if (cpList[i].size())
 899             {
 900                 file << cpList[i].size() << ' ';
 901
 902                 forAll(cpList[i], j)
 903                 {
 904                     file << cpList[i][j] << ' ';
 905                 }
 906
 907                 file << nl;
 908             }
 909         }
 910     }
 911     else
 912     {
 913         // List of points
 914         for (label i = 0; i < nPoints; i++)
 915         {
 916             file << 1 << ' ' << i << nl;
 917         }
 918     }
 919
 920     file << "CELL_TYPES " << nCells << endl;
 921
 922     if (cpList.size())
 923     {
 924         forAll(cpList, i)
 925         {
 926             switch (primitiveType)
 927             {
 928                 // Vertex
 929                 case 0:
 930                 {
 931                     file << "1" << nl;
 932                     break;
 933                 }
 934
 935                 // Edge
 936                 case 1:
 937                 {
 938                     file << "3" << nl;
 939                     break;
 940                 }
 941
 942                 // General Polygonal Face
 943                 case 2:
 944                 {
 945                     file << "7" << nl;
 946                     break;
 947                 }
 948
 949                 // Cell
 950                 case 3:
 951                 {
 952                     if (cpList[i].size() == 4)
 953                     {
 954                         // Tetrahedron
 955                         file << "10" << nl;
 956                     }
 957                     else
 958                     if (cpList[i].size() > 4)
 959                     {
 960                         // Polyhedron
 961                         file << "42" << nl;
 962                     }
 963
 964                     break;
 965                 }
 966
 967                 default:
 968                 {
 969                     FatalErrorIn("void meshOps::writeVTK() const")
 970                         << " Incorrect primitiveType: "
 971                         << primitiveType
 972                         << abort(FatalError);
 973                 }
 974             }
 975         }
 976     }
 977     else
 978     {
 979         // List of points
 980         for (label i = 0; i < nPoints; i++)
 981         {
 982             // Vertex
 983             file << '1' << nl;
 984         }
 985     }
 986
 987     label nCellFields = 0, nPointFields = 0;
 988
 989     if (reverseCellMap.size())
 990     {
 991         nCellFields++;
 992     }
 993
 994     if (scalField.size() == nCells && scalField.size())
 995     {
 996         nCellFields++;
 997     }
 998
 999     if (lablField.size() == nCells && lablField.size())
1000     {
1001         nCellFields++;
1002     }
1003
1004     if (vectField.size() == nCells && vectField.size())
1005     {
1006         nCellFields++;
1007     }
1008
1009     // Write out indices for visualization.
1010     if (nCellFields)
1011     {
1012         file << "CELL_DATA " << nCells << endl;
1013
1014         file << "FIELD CellFields " << nCellFields << endl;
1015
1016         if (reverseCellMap.size())
1017         {
1018             file << "CellIds 1 " << nCells << " int" << endl;
1019
1020             for (label i = 0; i < nCells; i++)
1021             {
1022                 file << reverseCellMap[i] << ' ';
1023             }
1024
1025             file << endl;
1026         }
1027
1028         if (scalField.size() == nCells)
1029         {
1030             file << "CellScalars 1 " << nCells << " double" << endl;
1031
1032             forAll(scalField, i)
1033             {
1034                 file << scalField[i] << ' ';
1035             }
1036
1037             file << endl;
1038         }
1039
1040         if (lablField.size() == nCells)
1041         {
1042             file << "CellLabels 1 " << nCells << " int" << endl;
1043
1044             forAll(lablField, i)
1045             {
1046                 file << lablField[i] << ' ';
1047             }
1048
1049             file << endl;
1050         }
1051
1052         if (vectField.size() == nCells)
1053         {
1054             file << "CellVectors 3 " << nCells << " double" << endl;
1055
1056             forAll(vectField, i)
1057             {
1058                 file << vectField[i].x() << ' '
1059                      << vectField[i].y() << ' '
1060                      << vectField[i].z() << ' ';
1061             }
1062
1063             file << endl;
1064         }
1065     }
1066
1067     if (reversePointMap.size())
1068     {
1069         nPointFields++;
1070     }
1071
1072     if (scalField.size() == nPoints && scalField.size())
1073     {
1074         nPointFields++;
1075     }
1076
1077     if (lablField.size() == nPoints && lablField.size())
1078     {
1079         nPointFields++;
1080     }
1081
1082     if (vectField.size() == nPoints && vectField.size())
1083     {
1084         nPointFields++;
1085     }
1086
1087     // Write out indices for visualization.
1088     if (nPointFields)
1089     {
1090         file << "POINT_DATA " << nPoints << endl;
1091
1092         file << "FIELD PointFields " << nPointFields << endl;
1093
1094         if (reversePointMap.size())
1095         {
1096             file << "PointIds 1 " << nPoints << " int" << endl;
1097
1098             for (label i = 0; i < nPoints; i++)
1099             {
1100                 file << reversePointMap[i] << ' ';
1101             }
1102
1103             file << endl;
1104         }
1105
1106         if (scalField.size() == nPoints)
1107         {
1108             file << "PointScalars 1 " << nPoints << " double" << endl;
1109
1110             forAll(scalField, i)
1111             {
1112                 file << scalField[i] << ' ';
1113             }
1114
1115             file << endl;
1116         }
1117
1118         if (lablField.size() == nPoints)
1119         {
1120             file << "PointLabels 1 " << nPoints << " int" << endl;
1121
1122             forAll(lablField, i)
1123             {
1124                 file << lablField[i] << ' ';
1125             }
1126
1127             file << endl;
1128         }
1129
1130         if (vectField.size() == nPoints)
1131         {
1132             file << "PointVectors 3 " << nPoints << " double" << endl;
1133
1134             forAll(vectField, i)
1135             {
1136                 file << vectField[i].x() << ' '
1137                      << vectField[i].y() << ' '
1138                      << vectField[i].z() << ' ';
1139             }
1140
1141             file << endl;
1142         }
1143     }
1144 }
1145
1146
1147 } // End namespace meshOps
1148
1149
1150 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1151
1152 } // End namespace Foam
1153
1154 // ************************************************************************* //