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