diff --git a/applications/utilities/postProcessing/graphics/PVReaders/PVFoamReader/vtkPVFoam/foamVtkAdaptors.H b/applications/utilities/postProcessing/graphics/PVReaders/PVFoamReader/vtkPVFoam/foamVtkAdaptors.H new file mode 100644 index 0000000000..89146bc889 --- /dev/null +++ b/applications/utilities/postProcessing/graphics/PVReaders/PVFoamReader/vtkPVFoam/foamVtkAdaptors.H @@ -0,0 +1,115 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2017 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#ifndef foamVtkAdaptors_H +#define foamVtkAdaptors_H + +// OpenFOAM includes +#include "labelList.H" + +// VTK includes +#include "vtkCellArray.h" +#include "vtkIdTypeArray.h" +#include "vtkSmartPointer.h" +#include "vtkUnsignedCharArray.h" +#include "vtkAOSDataArrayTemplate.h" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + //- Attach a smart pointer, or generate a non-null one. + template + inline vtkSmartPointer nonNullSmartPointer(T* ptr) + { + return vtkSmartPointer(ptr ? ptr : T::New()); + } + + + //- Helper to wrap vtkUnsignedCharArray as a UList + inline UList vtkUList + ( + vtkUnsignedCharArray* array, + const label size + ) + { + array->SetNumberOfComponents(1); + array->SetNumberOfTuples(size); + + UList list + ( + array->WritePointer(0, size), + size + ); + + return list; + } + + + //- Helper to wrap vtkIdTypeArray as a UList + inline UList vtkUList + ( + vtkIdTypeArray* array, + const label size + ) + { + array->SetNumberOfComponents(1); + array->SetNumberOfTuples(size); + + UList list + ( + array->WritePointer(0, size), + size + ); + + return list; + } + + + //- Special helper to wrap vtkCellArray as a UList + inline UList vtkUList + ( + vtkCellArray* cells, + const label nCells, + const label size + ) + { + cells->GetData()->SetNumberOfTuples(size); + + UList list + ( + cells->WritePointer(nCells, size), + size + ); + + return list; + } +} + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/utilities/postProcessing/graphics/PVReaders/PVFoamReader/vtkPVFoam/vtkPVFoamMeshVolume.C b/applications/utilities/postProcessing/graphics/PVReaders/PVFoamReader/vtkPVFoam/vtkPVFoamMeshVolume.C index 487687f689..48dc989ff3 100644 --- a/applications/utilities/postProcessing/graphics/PVReaders/PVFoamReader/vtkPVFoam/vtkPVFoamMeshVolume.C +++ b/applications/utilities/postProcessing/graphics/PVReaders/PVFoamReader/vtkPVFoam/vtkPVFoamMeshVolume.C @@ -28,13 +28,11 @@ License // OpenFOAM includes #include "fvMesh.H" -#include "cellModeller.H" +#include "foamVtuSizing.H" // VTK includes -#include "vtkCellArray.h" -#include "vtkIdTypeArray.h" #include "vtkUnstructuredGrid.h" -#include "vtkSmartPointer.h" +#include "foamVtkAdaptors.H" // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // @@ -44,383 +42,134 @@ vtkSmartPointer Foam::vtkPVFoam::volumeVTKMesh foamVtuData& vtuData ) { - const cellModel& tet = *(cellModeller::lookup("tet")); - const cellModel& pyr = *(cellModeller::lookup("pyr")); - const cellModel& prism = *(cellModeller::lookup("prism")); - const cellModel& wedge = *(cellModeller::lookup("wedge")); - const cellModel& tetWedge = *(cellModeller::lookup("tetWedge")); - const cellModel& hex = *(cellModeller::lookup("hex")); + if (debug) + { + Info<< " " << FUNCTION_NAME << endl; + printMemory(); + } + + foamVtuSizing sizing(mesh, !reader_->GetUseVTKPolyhedron()); vtkSmartPointer vtkmesh = vtkSmartPointer::New(); - if (debug) - { - Info<< " volumeVTKMesh" << endl; - printMemory(); - } + auto cellTypes = + nonNullSmartPointer(vtkmesh->GetCellTypesArray()); - const cellShapeList& cellShapes = mesh.cellShapes(); + auto cells = + nonNullSmartPointer(vtkmesh->GetCells()); - // Number of additional points needed by the decomposition of polyhedra - label nAddPoints = 0; + auto faces = + nonNullSmartPointer(vtkmesh->GetFaces()); - // Number of additional cells generated by the decomposition of polyhedra - label nAddCells = 0; + auto cellLocations = + nonNullSmartPointer(vtkmesh->GetCellLocationsArray()); - // face owner is needed to determine the face orientation - const labelList& owner = mesh.faceOwner(); + auto faceLocations = + nonNullSmartPointer(vtkmesh->GetFaceLocations()); - labelList& cellMap = vtuData.cellMap(); - labelList& addPointCellLabels = vtuData.additionalIds(); + UList cellTypesUL = + vtkUList + ( + cellTypes, + sizing.nFieldCells() + ); - // Scan for cells which need to be decomposed and count additional points - // and cells - if (!reader_->GetUseVTKPolyhedron()) - { - forAll(cellShapes, celli) - { - const cellModel& model = cellShapes[celli].model(); + UList cellsUL = + vtkUList + ( + cells, + sizing.nFieldCells(), + sizing.sizeInternal(foamVtuSizing::slotType::CELLS) + ); - if - ( - model != hex - && model != wedge - && model != prism - && model != pyr - && model != tet - && model != tetWedge - ) - { - const cell& cFaces = mesh.cells()[celli]; + UList cellLocationsUL = + vtkUList + ( + cellLocations, + sizing.sizeInternal(foamVtuSizing::slotType::CELLS_OFFSETS) + ); - forAll(cFaces, cFacei) - { - const face& f = mesh.faces()[cFaces[cFacei]]; + UList facesUL = + vtkUList + ( + faces, + sizing.sizeInternal(foamVtuSizing::slotType::FACES) + ); - label nQuads = 0; - label nTris = 0; - f.nTrianglesQuads(mesh.points(), nTris, nQuads); + UList faceLocationsUL = + vtkUList + ( + faceLocations, + sizing.sizeInternal(foamVtuSizing::slotType::FACES_OFFSETS) + ); - nAddCells += nQuads + nTris; - } - nAddCells--; - nAddPoints++; - } - } - } + sizing.populateInternal + ( + mesh, + cellTypesUL, + cellsUL, + cellLocationsUL, + facesUL, + faceLocationsUL, + static_cast(vtuData) + ); - // Set size of additional point addressing array - // (from added point to original cell) - addPointCellLabels.setSize(nAddPoints); - - // Set size of additional cells mapping array - // (from added cell to original cell) - - cellMap.setSize(mesh.nCells() + nAddCells); // Convert OpenFOAM mesh vertices to VTK - vtkSmartPointer vtkpoints = vtkSmartPointer::New(); + // - can only do this *after* populating the decompInfo with cell-ids + // for any additional points (ie, mesh cell-centres) + vtkSmartPointer vtkpoints = vtkmesh->GetPoints(); + if (!vtkpoints) + { + // No points previously, add an entry + vtkpoints = vtkSmartPointer::New(); + vtkmesh->SetPoints(vtkpoints); + } - vtkpoints->Allocate(mesh.nPoints() + nAddPoints); + vtkpoints->SetNumberOfPoints(sizing.nFieldPoints()); - const Foam::pointField& points = mesh.points(); + // Normal points + const pointField& points = mesh.points(); + const labelList& addPoints = vtuData.additionalIds(); + label pointi = 0; forAll(points, i) { - vtkpoints->InsertNextPoint(points[i].v_); + vtkpoints->SetPoint(pointi++, points[i].v_); } - - vtkmesh->Allocate(mesh.nCells() + nAddCells); - - // Set counters for additional points and additional cells - label addPointi = 0, addCelli = 0; - - // Create storage for points - needed for mapping from OpenFOAM to VTK - // data types - max 'order' = hex = 8 points - vtkIdType nodeIds[8]; - - // face-stream for a polyhedral cell - // [numFace0Pts, id1, id2, id3, numFace1Pts, id1, id2, id3, ...] - DynamicList faceStream(256); - - forAll(cellShapes, celli) + forAll(addPoints, i) { - const cellShape& cellShape = cellShapes[celli]; - const cellModel& cellModel = cellShape.model(); - - cellMap[addCelli++] = celli; - - if (cellModel == tet) - { - for (int j = 0; j < 4; j++) - { - nodeIds[j] = cellShape[j]; - } - vtkmesh->InsertNextCell - ( - VTK_TETRA, - 4, - nodeIds - ); - } - else if (cellModel == pyr) - { - for (int j = 0; j < 5; j++) - { - nodeIds[j] = cellShape[j]; - } - vtkmesh->InsertNextCell - ( - VTK_PYRAMID, - 5, - nodeIds - ); - } - else if (cellModel == prism) - { - // VTK has a different node order for VTK_WEDGE - // their triangles point outwards! - nodeIds[0] = cellShape[0]; - nodeIds[1] = cellShape[2]; - nodeIds[2] = cellShape[1]; - nodeIds[3] = cellShape[3]; - nodeIds[4] = cellShape[5]; - nodeIds[5] = cellShape[4]; - - vtkmesh->InsertNextCell - ( - VTK_WEDGE, - 6, - nodeIds - ); - } - else if (cellModel == tetWedge && !reader_->GetUseVTKPolyhedron()) - { - // Treat as squeezed prism (VTK_WEDGE) - - nodeIds[0] = cellShape[0]; - nodeIds[1] = cellShape[2]; - nodeIds[2] = cellShape[1]; - nodeIds[3] = cellShape[3]; - nodeIds[4] = cellShape[4]; - nodeIds[5] = cellShape[3]; - - vtkmesh->InsertNextCell - ( - VTK_WEDGE, - 6, - nodeIds - ); - } - else if (cellModel == wedge) - { - // Treat as squeezed hex - - nodeIds[0] = cellShape[0]; - nodeIds[1] = cellShape[1]; - nodeIds[2] = cellShape[2]; - nodeIds[3] = cellShape[2]; - nodeIds[4] = cellShape[3]; - nodeIds[5] = cellShape[4]; - nodeIds[6] = cellShape[5]; - nodeIds[7] = cellShape[6]; - - vtkmesh->InsertNextCell - ( - VTK_HEXAHEDRON, - 8, - nodeIds - ); - } - else if (cellModel == hex) - { - for (int j = 0; j < 8; j++) - { - nodeIds[j] = cellShape[j]; - } - vtkmesh->InsertNextCell - ( - VTK_HEXAHEDRON, - 8, - nodeIds - ); - } - else if (reader_->GetUseVTKPolyhedron()) - { - // Polyhedral cell - use VTK_POLYHEDRON - const labelList& cFaces = mesh.cells()[celli]; - - vtkIdType nFaces = cFaces.size(); - vtkIdType nLabels = nFaces; - - // count size for face stream - forAll(cFaces, cFacei) - { - const face& f = mesh.faces()[cFaces[cFacei]]; - nLabels += f.size(); - } - - // build face-stream - // [numFace0Pts, id1, id2, id3, numFace1Pts, id1, id2, id3, ...] - // point Ids are global - faceStream.clear(); - faceStream.reserve(nLabels + nFaces); - - forAll(cFaces, cFacei) - { - const face& f = mesh.faces()[cFaces[cFacei]]; - const bool isOwner = (owner[cFaces[cFacei]] == celli); - const label nFacePoints = f.size(); - - // number of labels for this face - faceStream.append(nFacePoints); - - if (isOwner) - { - forAll(f, fp) - { - faceStream.append(f[fp]); - } - } - else - { - // fairly immaterial if we reverse the list - // or use face::reverseFace() - forAllReverse(f, fp) - { - faceStream.append(f[fp]); - } - } - } - - vtkmesh->InsertNextCell(VTK_POLYHEDRON, nFaces, faceStream.data()); - } - else - { - // Polyhedral cell. Decompose into tets + prisms. - - // Mapping from additional point to cell - addPointCellLabels[addPointi] = celli; - - // The new vertex from the cell-centre - const label newVertexLabel = mesh.nPoints() + addPointi; - vtkpoints->InsertNextPoint(mesh.C()[celli].v_); - - // Whether to insert cell in place of original or not. - bool substituteCell = true; - - const labelList& cFaces = mesh.cells()[celli]; - forAll(cFaces, cFacei) - { - const face& f = mesh.faces()[cFaces[cFacei]]; - const bool isOwner = (owner[cFaces[cFacei]] == celli); - - // Number of triangles and quads in decomposition - label nTris = 0; - label nQuads = 0; - f.nTrianglesQuads(mesh.points(), nTris, nQuads); - - // Do actual decomposition into triFcs and quadFcs. - faceList triFcs(nTris); - faceList quadFcs(nQuads); - label trii = 0; - label quadi = 0; - f.trianglesQuads(mesh.points(), trii, quadi, triFcs, quadFcs); - - forAll(quadFcs, quadI) - { - if (substituteCell) - { - substituteCell = false; - } - else - { - cellMap[addCelli++] = celli; - } - - const face& quad = quadFcs[quadI]; - - // Ensure we have the correct orientation for the - // base of the primitive cell shape. - // If the cell is face owner, the orientation needs to be - // flipped. - // At the moment, VTK doesn't actually seem to care if - // negative cells are defined, but we'll do it anyhow - // (for safety). - if (isOwner) - { - nodeIds[0] = quad[3]; - nodeIds[1] = quad[2]; - nodeIds[2] = quad[1]; - nodeIds[3] = quad[0]; - } - else - { - nodeIds[0] = quad[0]; - nodeIds[1] = quad[1]; - nodeIds[2] = quad[2]; - nodeIds[3] = quad[3]; - } - nodeIds[4] = newVertexLabel; - vtkmesh->InsertNextCell - ( - VTK_PYRAMID, - 5, - nodeIds - ); - } - - forAll(triFcs, triI) - { - if (substituteCell) - { - substituteCell = false; - } - else - { - cellMap[addCelli++] = celli; - } - - const face& tri = triFcs[triI]; - - // See note above about the orientation. - if (isOwner) - { - nodeIds[0] = tri[2]; - nodeIds[1] = tri[1]; - nodeIds[2] = tri[0]; - } - else - { - nodeIds[0] = tri[0]; - nodeIds[1] = tri[1]; - nodeIds[2] = tri[2]; - } - nodeIds[3] = newVertexLabel; - - vtkmesh->InsertNextCell - ( - VTK_TETRA, - 4, - nodeIds - ); - } - } - - addPointi++; - } + vtkpoints->SetPoint(pointi++, mesh.C()[addPoints[i]].v_); } - vtkmesh->SetPoints(vtkpoints); + if (facesUL.size()) + { + vtkmesh->SetCells + ( + cellTypes, + cellLocations, + cells, + faceLocations, + faces + ); + } + else + { + vtkmesh->SetCells + ( + cellTypes, + cellLocations, + cells, + nullptr, + nullptr + ); + } if (debug) { - Info<<"nCells=" << mesh.nCells() <<" nPoints=" << mesh.nPoints() - <<" nAddCells=" << nAddCells <<" nAddPoints=" << nAddPoints - << nl - << " volumeVTKMesh" << endl; + Info<< " " << FUNCTION_NAME << endl; printMemory(); }