applications/utilities/postProcessing/graphics/PVFoamReader/vtkFoam/vtkFoamAddInternalMesh.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
   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 "vtkFoam.H"
  30 #include "fvMesh.H"
  31 #include "cellModeller.H"
  32
  33 #include "vtkUnstructuredGrid.h"
  34 #include "vtkCellArray.h"
  35
  36 #include "vtkFoamInsertNextPoint.H"
  37
  38 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  39
  40 void Foam::vtkFoam::addInternalMesh
  41 (
  42     const fvMesh& mesh,
  43     vtkUnstructuredGrid* vtkMesh
  44 )
  45 {
  46     SetName(vtkMesh, "Internal Mesh");
  47
  48     // Number of additional points needed by the decomposition of polyhedra
  49     label nAddPoints = 0;
  50
  51     // Number of additional cells generated by the decomposition of polyhedra
  52     label nAddCells = 0;
  53
  54     const cellModel& tet = *(cellModeller::lookup("tet"));
  55     const cellModel& pyr = *(cellModeller::lookup("pyr"));
  56     const cellModel& prism = *(cellModeller::lookup("prism"));
  57     const cellModel& wedge = *(cellModeller::lookup("wedge"));
  58     const cellModel& tetWedge = *(cellModeller::lookup("tetWedge"));
  59     const cellModel& hex = *(cellModeller::lookup("hex"));
  60
  61     // Scan for cells which need to be decomposed and count additional points
  62     // and cells
  63     if (debug)
  64     {
  65         Info<< "building cell-shapes" << endl;
  66     }
  67     const cellShapeList& cellShapes = mesh.cellShapes();
  68
  69     if (debug)
  70     {
  71         Info<< "scanning" << endl;
  72     }
  73
  74     forAll(cellShapes, cellI)
  75     {
  76         const cellModel& model = cellShapes[cellI].model();
  77
  78         if
  79         (
  80             model != hex
  81          && model != wedge
  82          && model != prism
  83          && model != pyr
  84          && model != tet
  85          && model != tetWedge
  86         )
  87         {
  88             const cell& cFaces = mesh.cells()[cellI];
  89
  90             forAll(cFaces, cFaceI)
  91             {
  92                 const face& f = mesh.faces()[cFaces[cFaceI]];
  93
  94                 label nFacePoints = f.size();
  95
  96                 label nQuads = (nFacePoints - 2)/2;
  97                 label nTris = (nFacePoints - 2)%2;
  98                 nAddCells += nQuads + nTris;
  99             }
 100
 101             nAddCells--;
 102             nAddPoints++;
 103         }
 104     }
 105
 106     // Set size of additional point addressing array
 107     // (from added point to original cell)
 108     addPointCellLabels_.setSize(nAddPoints);
 109
 110     // Set size of additional cells mapping array
 111     // (from added cell to original cell)
 112     superCells_.setSize(mesh.nCells() + nAddCells);
 113
 114     if (debug)
 115     {
 116         Info<< "converting points" << endl;
 117     }
 118
 119     // Convert Foam mesh vertices to VTK
 120     vtkPoints *vtkpoints = vtkPoints::New();
 121     vtkpoints->Allocate(mesh.nPoints() + nAddPoints);
 122
 123     const Foam::pointField& points = mesh.points();
 124
 125     forAll(points, i)
 126     {
 127         vtkFoamInsertNextPoint(vtkpoints, points[i]);
 128     }
 129
 130     if (debug)
 131     {
 132         Info<< "converting cells" << endl;
 133     }
 134
 135     vtkMesh->Allocate(mesh.nCells() + nAddCells);
 136
 137     // Set counters for additional points and additional cells
 138     label api = 0, aci = 0;
 139
 140     forAll(cellShapes, celli)
 141     {
 142         const cellShape& cellShape = cellShapes[celli];
 143         const cellModel& cellModel = cellShape.model();
 144
 145         superCells_[aci++] = celli;
 146
 147         if (cellModel == tet)
 148         {
 149             vtkMesh->InsertNextCell
 150             (
 151                 VTK_TETRA,
 152                 4,
 153                 const_cast<int*>(cellShape.begin())
 154             );
 155         }
 156         else if (cellModel == pyr)
 157         {
 158             vtkMesh->InsertNextCell
 159             (
 160                 VTK_PYRAMID,
 161                 5,
 162                 const_cast<int*>(cellShape.begin())
 163             );
 164         }
 165         else if (cellModel == prism)
 166         {
 167             vtkMesh->InsertNextCell
 168             (
 169                 VTK_WEDGE,
 170                 6,
 171                 const_cast<int*>(cellShape.begin())
 172             );
 173         }
 174         else if (cellModel == tetWedge)
 175         {
 176             // Treat as squeezed prism
 177
 178             int vtkVerts[6];
 179             vtkVerts[0] = cellShape[0];
 180             vtkVerts[1] = cellShape[2];
 181             vtkVerts[2] = cellShape[1];
 182             vtkVerts[3] = cellShape[3];
 183             vtkVerts[4] = cellShape[4];
 184             vtkVerts[5] = cellShape[4];
 185
 186             vtkMesh->InsertNextCell(VTK_WEDGE, 6, vtkVerts);
 187         }
 188         else if (cellModel == wedge)
 189         {
 190             // Treat as squeezed hex
 191
 192             int vtkVerts[8];
 193             vtkVerts[0] = cellShape[0];
 194             vtkVerts[1] = cellShape[1];
 195             vtkVerts[2] = cellShape[2];
 196             vtkVerts[3] = cellShape[2];
 197             vtkVerts[4] = cellShape[3];
 198             vtkVerts[5] = cellShape[4];
 199             vtkVerts[6] = cellShape[5];
 200             vtkVerts[7] = cellShape[6];
 201
 202             vtkMesh->InsertNextCell(VTK_HEXAHEDRON, 8, vtkVerts);
 203         }
 204         else if (cellModel == hex)
 205         {
 206             vtkMesh->InsertNextCell
 207             (
 208                 VTK_HEXAHEDRON,
 209                 8,
 210                 const_cast<int*>(cellShape.begin())
 211             );
 212         }
 213         else
 214         {
 215             // Polyhedral cell. Decompose into tets + prisms.
 216
 217             // Mapping from additional point to cell
 218             addPointCellLabels_[api] = celli;
 219
 220             // Insert the new vertex from the cell-centre
 221             label newVertexLabel = mesh.nPoints() + api;
 222             vtkFoamInsertNextPoint(vtkpoints, mesh.C()[celli]);
 223
 224             // Whether to insert cell in place of original or not.
 225             bool substituteCell = true;
 226
 227             const labelList& cFaces = mesh.cells()[celli];
 228
 229             forAll(cFaces, cFaceI)
 230             {
 231                 const face& f = mesh.faces()[cFaces[cFaceI]];
 232
 233                 label nFacePoints = f.size();
 234
 235                 label nQuads = (nFacePoints - 2)/2;
 236                 label nTris = (nFacePoints - 2)%2;
 237
 238                 label qpi = 0;
 239
 240                 for (label quadi=0; quadi<nQuads; quadi++)
 241                 {
 242                     label thisCellI = -1;
 243
 244                     if (substituteCell)
 245                     {
 246                         thisCellI = celli;
 247                         substituteCell = false;
 248                     }
 249                     else
 250                     {
 251                         thisCellI = mesh.nCells() + aci;
 252                         superCells_[aci++] = celli;
 253                     }
 254
 255                     int addVtkVerts[5];
 256                     addVtkVerts[0] = f[0];
 257                     addVtkVerts[1] = f[qpi + 1];
 258                     addVtkVerts[2] = f[qpi + 2];
 259                     addVtkVerts[3] = f[qpi + 3];
 260                     addVtkVerts[4] = newVertexLabel;
 261                     vtkMesh->InsertNextCell(VTK_PYRAMID, 5, addVtkVerts);
 262
 263                     qpi += 2;
 264                 }
 265
 266                 if (nTris)
 267                 {
 268                     label thisCellI = -1;
 269
 270                     if (substituteCell)
 271                     {
 272                         thisCellI = celli;
 273                         substituteCell = false;
 274                     }
 275                     else
 276                     {
 277                         thisCellI = mesh.nCells() + aci;
 278                         superCells_[aci++] = celli;
 279                     }
 280
 281                     int addVtkVerts[4];
 282                     addVtkVerts[0] = f[0];
 283                     addVtkVerts[1] = f[qpi + 1];
 284                     addVtkVerts[2] = f[qpi + 2];
 285                     addVtkVerts[3] = newVertexLabel;
 286                     vtkMesh->InsertNextCell(VTK_TETRA, 4, addVtkVerts);
 287                 }
 288             }
 289
 290             api++;
 291         }
 292     }
 293
 294     vtkMesh->SetPoints(vtkpoints);
 295     vtkpoints->Delete();
 296 }
 297
 298
 299 // ************************************************************************* //