applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamMeshVolume.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   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
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your 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, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "vtkPV3Foam.H"
  27 #include "vtkPV3FoamReader.h"
  28
  29 // OpenFOAM includes
  30 #include "fvMesh.H"
  31 #include "cellModeller.H"
  32 #include "vtkOpenFOAMPoints.H"
  33 #include "Swap.H"
  34 #include "longLong.H"
  35
  36 // VTK includes
  37 #include "vtkCellArray.h"
  38 #include "vtkIdTypeArray.h"
  39 #include "vtkUnstructuredGrid.h"
  40
  41 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  42
  43 vtkUnstructuredGrid* Foam::vtkPV3Foam::volumeVTKMesh
  44 (
  45     const fvMesh& mesh,
  46     polyDecomp& decompInfo
  47 )
  48 {
  49     const cellModel& tet = *(cellModeller::lookup("tet"));
  50     const cellModel& pyr = *(cellModeller::lookup("pyr"));
  51     const cellModel& prism = *(cellModeller::lookup("prism"));
  52     const cellModel& wedge = *(cellModeller::lookup("wedge"));
  53     const cellModel& tetWedge = *(cellModeller::lookup("tetWedge"));
  54     const cellModel& hex = *(cellModeller::lookup("hex"));
  55
  56     vtkUnstructuredGrid* vtkmesh = vtkUnstructuredGrid::New();
  57
  58     if (debug)
  59     {
  60         Info<< "<beg> Foam::vtkPV3Foam::volumeVTKMesh" << endl;
  61         printMemory();
  62     }
  63
  64     const cellShapeList& cellShapes = mesh.cellShapes();
  65
  66     // Number of additional points needed by the decomposition of polyhedra
  67     label nAddPoints = 0;
  68
  69     // Number of additional cells generated by the decomposition of polyhedra
  70     label nAddCells = 0;
  71
  72     // face owner is needed to determine the face orientation
  73     const labelList& owner = mesh.faceOwner();
  74
  75     labelList& superCells = decompInfo.superCells();
  76     labelList& addPointCellLabels = decompInfo.addPointCellLabels();
  77
  78     // Scan for cells which need to be decomposed and count additional points
  79     // and cells
  80     if (!reader_->GetUseVTKPolyhedron())
  81     {
  82         if (debug)
  83         {
  84             Info<< "... scanning for polyhedra" << endl;
  85         }
  86
  87         forAll(cellShapes, cellI)
  88         {
  89             const cellModel& model = cellShapes[cellI].model();
  90
  91             if
  92             (
  93                 model != hex
  94              && model != wedge
  95              && model != prism
  96              && model != pyr
  97              && model != tet
  98              && model != tetWedge
  99             )
 100             {
 101                 const cell& cFaces = mesh.cells()[cellI];
 102
 103                 forAll(cFaces, cFaceI)
 104                 {
 105                     const face& f = mesh.faces()[cFaces[cFaceI]];
 106
 107                     label nQuads = 0;
 108                     label nTris = 0;
 109                     f.nTrianglesQuads(mesh.points(), nTris, nQuads);
 110
 111                     nAddCells += nQuads + nTris;
 112                 }
 113
 114                 nAddCells--;
 115                 nAddPoints++;
 116             }
 117         }
 118     }
 119
 120     // Set size of additional point addressing array
 121     // (from added point to original cell)
 122     addPointCellLabels.setSize(nAddPoints);
 123
 124     // Set size of additional cells mapping array
 125     // (from added cell to original cell)
 126
 127     if (debug)
 128     {
 129         Info<<" mesh nCells     = " << mesh.nCells() << nl
 130             <<"      nPoints    = " << mesh.nPoints() << nl
 131             <<"      nAddCells  = " << nAddCells << nl
 132             <<"      nAddPoints = " << nAddPoints << endl;
 133     }
 134
 135     superCells.setSize(mesh.nCells() + nAddCells);
 136
 137     if (debug)
 138     {
 139         Info<< "... converting points" << endl;
 140     }
 141
 142     // Convert OpenFOAM mesh vertices to VTK
 143     vtkPoints* vtkpoints = vtkPoints::New();
 144     vtkpoints->Allocate(mesh.nPoints() + nAddPoints);
 145
 146     const Foam::pointField& points = mesh.points();
 147
 148     forAll(points, i)
 149     {
 150         vtkInsertNextOpenFOAMPoint(vtkpoints, points[i]);
 151     }
 152
 153
 154     if (debug)
 155     {
 156         Info<< "... converting cells" << endl;
 157     }
 158
 159     vtkmesh->Allocate(mesh.nCells() + nAddCells);
 160
 161     // Set counters for additional points and additional cells
 162     label addPointI = 0, addCellI = 0;
 163
 164     // Create storage for points - needed for mapping from OpenFOAM to VTK
 165     // data types - max 'order' = hex = 8 points
 166     vtkIdType nodeIds[8];
 167
 168     // face-stream for a polyhedral cell
 169     // [numFace0Pts, id1, id2, id3, numFace1Pts, id1, id2, id3, ...]
 170     DynamicList<vtkIdType> faceStream(256);
 171
 172     forAll(cellShapes, cellI)
 173     {
 174         const cellShape& cellShape = cellShapes[cellI];
 175         const cellModel& cellModel = cellShape.model();
 176
 177         superCells[addCellI++] = cellI;
 178
 179         if (cellModel == tet)
 180         {
 181             for (int j = 0; j < 4; j++)
 182             {
 183                 nodeIds[j] = cellShape[j];
 184             }
 185             vtkmesh->InsertNextCell
 186             (
 187                 VTK_TETRA,
 188                 4,
 189                 nodeIds
 190             );
 191         }
 192         else if (cellModel == pyr)
 193         {
 194             for (int j = 0; j < 5; j++)
 195             {
 196                 nodeIds[j] = cellShape[j];
 197             }
 198             vtkmesh->InsertNextCell
 199             (
 200                 VTK_PYRAMID,
 201                 5,
 202                 nodeIds
 203             );
 204         }
 205         else if (cellModel == prism)
 206         {
 207             // VTK has a different node order for VTK_WEDGE
 208             // their triangles point outwards!
 209             nodeIds[0] = cellShape[0];
 210             nodeIds[1] = cellShape[2];
 211             nodeIds[2] = cellShape[1];
 212             nodeIds[3] = cellShape[3];
 213             nodeIds[4] = cellShape[5];
 214             nodeIds[5] = cellShape[4];
 215
 216             vtkmesh->InsertNextCell
 217             (
 218                 VTK_WEDGE,
 219                 6,
 220                 nodeIds
 221             );
 222         }
 223         else if (cellModel == tetWedge)
 224         {
 225             // Treat as squeezed prism
 226
 227             nodeIds[0] = cellShape[0];
 228             nodeIds[1] = cellShape[2];
 229             nodeIds[2] = cellShape[1];
 230             nodeIds[3] = cellShape[3];
 231             nodeIds[4] = cellShape[4];
 232             nodeIds[5] = cellShape[4];
 233
 234             vtkmesh->InsertNextCell
 235             (
 236                 VTK_WEDGE,
 237                 6,
 238                 nodeIds
 239             );
 240         }
 241         else if (cellModel == wedge)
 242         {
 243             // Treat as squeezed hex
 244
 245             nodeIds[0] = cellShape[0];
 246             nodeIds[1] = cellShape[1];
 247             nodeIds[2] = cellShape[2];
 248             nodeIds[3] = cellShape[2];
 249             nodeIds[4] = cellShape[3];
 250             nodeIds[5] = cellShape[4];
 251             nodeIds[6] = cellShape[5];
 252             nodeIds[7] = cellShape[6];
 253
 254             vtkmesh->InsertNextCell
 255             (
 256                 VTK_HEXAHEDRON,
 257                 8,
 258                 nodeIds
 259             );
 260         }
 261         else if (cellModel == hex)
 262         {
 263             for (int j = 0; j < 8; j++)
 264             {
 265                 nodeIds[j] = cellShape[j];
 266             }
 267             vtkmesh->InsertNextCell
 268             (
 269                 VTK_HEXAHEDRON,
 270                 8,
 271                 nodeIds
 272             );
 273         }
 274         else if (reader_->GetUseVTKPolyhedron())
 275         {
 276             // Polyhedral cell - use VTK_POLYHEDRON
 277             const labelList& cFaces = mesh.cells()[cellI];
 278
 279 #ifdef HAS_VTK_POLYHEDRON
 280             vtkIdType nFaces = cFaces.size();
 281             vtkIdType nLabels = nFaces;
 282
 283             // count size for face stream
 284             forAll(cFaces, cFaceI)
 285             {
 286                 const face& f = mesh.faces()[cFaces[cFaceI]];
 287                 nLabels += f.size();
 288             }
 289
 290             // build face-stream
 291             // [numFace0Pts, id1, id2, id3, numFace1Pts, id1, id2, id3, ...]
 292             // point Ids are global
 293             faceStream.clear();
 294             faceStream.reserve(nLabels + nFaces);
 295
 296             forAll(cFaces, cFaceI)
 297             {
 298                 const face& f = mesh.faces()[cFaces[cFaceI]];
 299                 const bool isOwner = (owner[cFaces[cFaceI]] == cellI);
 300                 const label nFacePoints = f.size();
 301
 302                 // number of labels for this face
 303                 faceStream.append(nFacePoints);
 304
 305                 if (isOwner)
 306                 {
 307                     forAll(f, fp)
 308                     {
 309                         faceStream.append(f[fp]);
 310                     }
 311                 }
 312                 else
 313                 {
 314                     // fairly immaterial if we reverse the list
 315                     // or use face::reverseFace()
 316                     forAllReverse(f, fp)
 317                     {
 318                         faceStream.append(f[fp]);
 319                     }
 320                 }
 321             }
 322
 323             vtkmesh->InsertNextCell(VTK_POLYHEDRON, nFaces, faceStream.data());
 324 #else
 325             // this is a horrible substitute
 326             // but avoids crashes when there is no vtkPolyhedron support
 327
 328             // establish unique node ids used
 329             HashSet<vtkIdType, Hash<label> > hashUniqId(2*256);
 330
 331             forAll(cFaces, cFaceI)
 332             {
 333                 const face& f = mesh.faces()[cFaces[cFaceI]];
 334
 335                 forAll(f, fp)
 336                 {
 337                     hashUniqId.insert(f[fp]);
 338                 }
 339             }
 340
 341             // use face stream to store unique node ids:
 342             faceStream = hashUniqId.sortedToc();
 343
 344             vtkmesh->InsertNextCell
 345             (
 346                 VTK_CONVEX_POINT_SET,
 347                 vtkIdType(faceStream.size()),
 348                 faceStream.data()
 349             );
 350 #endif
 351         }
 352         else
 353         {
 354             // Polyhedral cell. Decompose into tets + prisms.
 355
 356             // Mapping from additional point to cell
 357             addPointCellLabels[addPointI] = cellI;
 358
 359             // The new vertex from the cell-centre
 360             const label newVertexLabel = mesh.nPoints() + addPointI;
 361             vtkInsertNextOpenFOAMPoint(vtkpoints, mesh.C()[cellI]);
 362
 363             // Whether to insert cell in place of original or not.
 364             bool substituteCell = true;
 365
 366             const labelList& cFaces = mesh.cells()[cellI];
 367             forAll(cFaces, cFaceI)
 368             {
 369                 const face& f = mesh.faces()[cFaces[cFaceI]];
 370                 const bool isOwner = (owner[cFaces[cFaceI]] == cellI);
 371
 372                 // Number of triangles and quads in decomposition
 373                 label nTris = 0;
 374                 label nQuads = 0;
 375                 f.nTrianglesQuads(mesh.points(), nTris, nQuads);
 376
 377                 // Do actual decomposition into triFcs and quadFcs.
 378                 faceList triFcs(nTris);
 379                 faceList quadFcs(nQuads);
 380                 label trii = 0;
 381                 label quadi = 0;
 382                 f.trianglesQuads(mesh.points(), trii, quadi, triFcs, quadFcs);
 383
 384                 forAll(quadFcs, quadI)
 385                 {
 386                     if (substituteCell)
 387                     {
 388                         substituteCell = false;
 389                     }
 390                     else
 391                     {
 392                         superCells[addCellI++] = cellI;
 393                     }
 394
 395                     const face& quad = quadFcs[quadI];
 396
 397                     // Ensure we have the correct orientation for the
 398                     // base of the primitive cell shape.
 399                     // If the cell is face owner, the orientation needs to be
 400                     // flipped.
 401                     // At the moment, VTK doesn't actually seem to care if
 402                     // negative cells are defined, but we'll do it anyhow
 403                     // (for safety).
 404                     if (isOwner)
 405                     {
 406                         nodeIds[0] = quad[3];
 407                         nodeIds[1] = quad[2];
 408                         nodeIds[2] = quad[1];
 409                         nodeIds[3] = quad[0];
 410                     }
 411                     else
 412                     {
 413                         nodeIds[0] = quad[0];
 414                         nodeIds[1] = quad[1];
 415                         nodeIds[2] = quad[2];
 416                         nodeIds[3] = quad[3];
 417                     }
 418                     nodeIds[4] = newVertexLabel;
 419                     vtkmesh->InsertNextCell
 420                     (
 421                         VTK_PYRAMID,
 422                         5,
 423                         nodeIds
 424                     );
 425                 }
 426
 427                 forAll(triFcs, triI)
 428                 {
 429                     if (substituteCell)
 430                     {
 431                         substituteCell = false;
 432                     }
 433                     else
 434                     {
 435                         superCells[addCellI++] = cellI;
 436                     }
 437
 438                     const face& tri = triFcs[triI];
 439
 440                     // See note above about the orientation.
 441                     if (isOwner)
 442                     {
 443                         nodeIds[0] = tri[2];
 444                         nodeIds[1] = tri[1];
 445                         nodeIds[2] = tri[0];
 446                     }
 447                     else
 448                     {
 449                         nodeIds[0] = tri[0];
 450                         nodeIds[1] = tri[1];
 451                         nodeIds[2] = tri[2];
 452                     }
 453                     nodeIds[3] = newVertexLabel;
 454
 455                     vtkmesh->InsertNextCell
 456                     (
 457                         VTK_TETRA,
 458                         4,
 459                         nodeIds
 460                     );
 461                 }
 462             }
 463
 464             addPointI++;
 465         }
 466     }
 467
 468     vtkmesh->SetPoints(vtkpoints);
 469     vtkpoints->Delete();
 470
 471     if (debug)
 472     {
 473         Info<< "<end> Foam::vtkPV3Foam::volumeVTKMesh" << endl;
 474         printMemory();
 475     }
 476
 477     return vtkmesh;
 478 }
 479
 480
 481 // ************************************************************************* //