applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamMeshVolume.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 Description
  26
  27 \*---------------------------------------------------------------------------*/
  28
  29 #include "vtkPV3Foam.H"
  30
  31 // Foam includes
  32 #include "fvMesh.H"
  33 #include "cellModeller.H"
  34 #include "vtkPV3FoamPoints.H"
  35
  36 // VTK includes
  37 #include "vtkCellArray.h"
  38 #include "vtkUnstructuredGrid.h"
  39
  40 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  41
  42 vtkUnstructuredGrid* Foam::vtkPV3Foam::volumeVTKMesh
  43 (
  44     const fvMesh& mesh,
  45     polyDecomp& decompInfo
  46 )
  47 {
  48     vtkUnstructuredGrid* vtkmesh = vtkUnstructuredGrid::New();
  49
  50     if (debug)
  51     {
  52         Info<< "<beg> Foam::vtkPV3Foam::volumeVTKMesh" << endl;
  53         printMemory();
  54     }
  55
  56     // Number of additional points needed by the decomposition of polyhedra
  57     label nAddPoints = 0;
  58
  59     // Number of additional cells generated by the decomposition of polyhedra
  60     label nAddCells = 0;
  61
  62     labelList& superCells = decompInfo.superCells();
  63     labelList& addPointCellLabels = decompInfo.addPointCellLabels();
  64
  65     const cellModel& tet = *(cellModeller::lookup("tet"));
  66     const cellModel& pyr = *(cellModeller::lookup("pyr"));
  67     const cellModel& prism = *(cellModeller::lookup("prism"));
  68     const cellModel& wedge = *(cellModeller::lookup("wedge"));
  69     const cellModel& tetWedge = *(cellModeller::lookup("tetWedge"));
  70     const cellModel& hex = *(cellModeller::lookup("hex"));
  71
  72     // Scan for cells which need to be decomposed and count additional points
  73     // and cells
  74     if (debug)
  75     {
  76         Info<< "... building cell-shapes" << endl;
  77     }
  78     const cellShapeList& cellShapes = mesh.cellShapes();
  79
  80     if (debug)
  81     {
  82         Info<< "... scanning" << endl;
  83     }
  84     forAll(cellShapes, cellI)
  85     {
  86         const cellModel& model = cellShapes[cellI].model();
  87
  88         if
  89         (
  90             model != hex
  91          && model != wedge
  92          && model != prism
  93          && model != pyr
  94          && model != tet
  95          && model != tetWedge
  96         )
  97         {
  98             const cell& cFaces = mesh.cells()[cellI];
  99
 100             forAll(cFaces, cFaceI)
 101             {
 102                 const face& f = mesh.faces()[cFaces[cFaceI]];
 103
 104                 label nFacePoints = f.size();
 105
 106                 label nQuads = (nFacePoints - 2)/2;
 107                 label nTris = (nFacePoints - 2)%2;
 108                 nAddCells += nQuads + nTris;
 109             }
 110
 111             nAddCells--;
 112             nAddPoints++;
 113         }
 114     }
 115
 116     // Set size of additional point addressing array
 117     // (from added point to original cell)
 118     addPointCellLabels.setSize(nAddPoints);
 119
 120     // Set size of additional cells mapping array
 121     // (from added cell to original cell)
 122
 123     if (debug)
 124     {
 125         Info<<" mesh nCells     = " << mesh.nCells() << nl
 126             <<"      nPoints    = " << mesh.nPoints() << nl
 127             <<"      nAddCells  = " << nAddCells << nl
 128             <<"      nAddPoints = " << nAddPoints << endl;
 129     }
 130
 131     superCells.setSize(mesh.nCells() + nAddCells);
 132
 133     if (debug)
 134     {
 135         Info<< "... converting points" << endl;
 136     }
 137
 138     // Convert Foam mesh vertices to VTK
 139     vtkPoints *vtkpoints = vtkPoints::New();
 140     vtkpoints->Allocate( mesh.nPoints() + nAddPoints );
 141
 142     const Foam::pointField& points = mesh.points();
 143
 144     forAll(points, i)
 145     {
 146         vtkPV3FoamInsertNextPoint(vtkpoints, points[i]);
 147     }
 148
 149
 150     if (debug)
 151     {
 152         Info<< "... converting cells" << endl;
 153     }
 154
 155     vtkmesh->Allocate( mesh.nCells() + nAddCells );
 156
 157     // Set counters for additional points and additional cells
 158     label addPointI = 0, addCellI = 0;
 159
 160     // Create storage for points - needed for mapping from Foam to VTK
 161     // data types - max 'order' = hex = 8 points
 162     vtkIdType nodeIds[8];
 163
 164     forAll(cellShapes, cellI)
 165     {
 166         const cellShape& cellShape = cellShapes[cellI];
 167         const cellModel& cellModel = cellShape.model();
 168
 169         superCells[addCellI++] = cellI;
 170
 171         if (cellModel == tet)
 172         {
 173             for (int j = 0; j < 4; j++)
 174             {
 175                 nodeIds[j] = cellShape[j];
 176             }
 177             vtkmesh->InsertNextCell
 178             (
 179                 VTK_TETRA,
 180                 4,
 181                 nodeIds
 182             );
 183         }
 184         else if (cellModel == pyr)
 185         {
 186             for (int j = 0; j < 5; j++)
 187             {
 188                 nodeIds[j] = cellShape[j];
 189             }
 190             vtkmesh->InsertNextCell
 191             (
 192                 VTK_PYRAMID,
 193                 5,
 194                 nodeIds
 195             );
 196         }
 197         else if (cellModel == prism)
 198         {
 199             for (int j = 0; j < 6; j++)
 200             {
 201                 nodeIds[j] = cellShape[j];
 202             }
 203             vtkmesh->InsertNextCell
 204             (
 205                 VTK_WEDGE,
 206                 6,
 207                 nodeIds
 208             );
 209         }
 210         else if (cellModel == tetWedge)
 211         {
 212             // Treat as squeezed prism
 213
 214             nodeIds[0] = cellShape[0];
 215             nodeIds[1] = cellShape[2];
 216             nodeIds[2] = cellShape[1];
 217             nodeIds[3] = cellShape[3];
 218             nodeIds[4] = cellShape[4];
 219             nodeIds[5] = cellShape[4];
 220
 221             vtkmesh->InsertNextCell
 222             (
 223                 VTK_WEDGE,
 224                 6,
 225                 nodeIds
 226             );
 227         }
 228         else if (cellModel == wedge)
 229         {
 230             // Treat as squeezed hex
 231
 232             nodeIds[0] = cellShape[0];
 233             nodeIds[1] = cellShape[1];
 234             nodeIds[2] = cellShape[2];
 235             nodeIds[3] = cellShape[2];
 236             nodeIds[4] = cellShape[3];
 237             nodeIds[5] = cellShape[4];
 238             nodeIds[6] = cellShape[5];
 239             nodeIds[7] = cellShape[6];
 240
 241             vtkmesh->InsertNextCell
 242             (
 243                 VTK_HEXAHEDRON,
 244                 8,
 245                 nodeIds
 246             );
 247         }
 248         else if (cellModel == hex)
 249         {
 250             for (int j = 0; j < 8; j++)
 251             {
 252                 nodeIds[j] = cellShape[j];
 253             }
 254             vtkmesh->InsertNextCell
 255             (
 256                 VTK_HEXAHEDRON,
 257                 8,
 258                 nodeIds
 259             );
 260         }
 261         else
 262         {
 263             // Polyhedral cell. Decompose into tets + prisms.
 264
 265             // Mapping from additional point to cell
 266             addPointCellLabels[addPointI] = cellI;
 267
 268             // Insert the new vertex from the cell-centre
 269             label newVertexLabel = mesh.nPoints() + addPointI;
 270             vtkPV3FoamInsertNextPoint(vtkpoints, mesh.C()[cellI]);
 271
 272             // Whether to insert cell in place of original or not.
 273             bool substituteCell = true;
 274
 275             const labelList& cFaces = mesh.cells()[cellI];
 276
 277             forAll(cFaces, cFaceI)
 278             {
 279                 const face& f = mesh.faces()[cFaces[cFaceI]];
 280
 281                 label nFacePoints = f.size();
 282
 283                 label nQuads = (nFacePoints - 2)/2;
 284                 label nTris = (nFacePoints - 2)%2;
 285
 286                 label qpi = 0;
 287
 288                 for (label quadi=0; quadi<nQuads; quadi++)
 289                 {
 290                     label thisCellI = -1;
 291
 292                     if (substituteCell)
 293                     {
 294                         thisCellI = cellI;
 295                         substituteCell = false;
 296                     }
 297                     else
 298                     {
 299                         thisCellI = mesh.nCells() + addCellI;
 300                         superCells[addCellI++] = cellI;
 301                     }
 302
 303                     nodeIds[0] = f[0];
 304                     nodeIds[1] = f[qpi + 1];
 305                     nodeIds[2] = f[qpi + 2];
 306                     nodeIds[3] = f[qpi + 3];
 307                     nodeIds[4] = newVertexLabel;
 308                     vtkmesh->InsertNextCell
 309                     (
 310                         VTK_PYRAMID,
 311                         5,
 312                         nodeIds
 313                     );
 314
 315                     qpi += 2;
 316                 }
 317
 318                 if (nTris)
 319                 {
 320                     label thisCellI = -1;
 321
 322                     if (substituteCell)
 323                     {
 324                         thisCellI = cellI;
 325                         substituteCell = false;
 326                     }
 327                     else
 328                     {
 329                         thisCellI = mesh.nCells() + addCellI;
 330                         superCells[addCellI++] = cellI;
 331                     }
 332
 333                     nodeIds[0] = f[0];
 334                     nodeIds[1] = f[qpi + 1];
 335                     nodeIds[2] = f[qpi + 2];
 336                     nodeIds[3] = newVertexLabel;
 337                     vtkmesh->InsertNextCell
 338                     (
 339                         VTK_TETRA,
 340                         4,
 341                         nodeIds
 342                     );
 343                 }
 344             }
 345
 346             addPointI++;
 347         }
 348     }
 349
 350     vtkmesh->SetPoints(vtkpoints);
 351     vtkpoints->Delete();
 352
 353     if (debug)
 354     {
 355         Info<< "<end> Foam::vtkPV3Foam::volumeVTKMesh" << endl;
 356         printMemory();
 357     }
 358
 359     return vtkmesh;
 360 }
 361
 362
 363 // ************************************************************************* //