Commit 8a66ec99 authored by Mark OLESEN's avatar Mark OLESEN
Browse files

ENH: replace paraview vtk conversion with library utilities

- this allows filling in the VTK structures without intermediate data
  and without sequencial insertion. Should be faster and smaller
  than the previous cell-wise insertion methods.
  Most importantly, it improves code reuse.
parent 1f18a0f9
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#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<class T>
inline vtkSmartPointer<T> nonNullSmartPointer(T* ptr)
{
return vtkSmartPointer<T>(ptr ? ptr : T::New());
}
//- Helper to wrap vtkUnsignedCharArray as a UList
inline UList<uint8_t> vtkUList
(
vtkUnsignedCharArray* array,
const label size
)
{
array->SetNumberOfComponents(1);
array->SetNumberOfTuples(size);
UList<uint8_t> list
(
array->WritePointer(0, size),
size
);
return list;
}
//- Helper to wrap vtkIdTypeArray as a UList
inline UList<vtkIdType> vtkUList
(
vtkIdTypeArray* array,
const label size
)
{
array->SetNumberOfComponents(1);
array->SetNumberOfTuples(size);
UList<vtkIdType> list
(
array->WritePointer(0, size),
size
);
return list;
}
//- Special helper to wrap vtkCellArray as a UList
inline UList<vtkIdType> vtkUList
(
vtkCellArray* cells,
const label nCells,
const label size
)
{
cells->GetData()->SetNumberOfTuples(size);
UList<vtkIdType> list
(
cells->WritePointer(nCells, size),
size
);
return list;
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //
......@@ -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<vtkUnstructuredGrid> 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"));
vtkSmartPointer<vtkUnstructuredGrid> vtkmesh =
vtkSmartPointer<vtkUnstructuredGrid>::New();
if (debug)
{
Info<< "<beg> volumeVTKMesh" << endl;
Info<< "<beg> " << FUNCTION_NAME << endl;
printMemory();
}
const cellShapeList& cellShapes = mesh.cellShapes();
// Number of additional points needed by the decomposition of polyhedra
label nAddPoints = 0;
foamVtuSizing sizing(mesh, !reader_->GetUseVTKPolyhedron());
// Number of additional cells generated by the decomposition of polyhedra
label nAddCells = 0;
vtkSmartPointer<vtkUnstructuredGrid> vtkmesh =
vtkSmartPointer<vtkUnstructuredGrid>::New();
// face owner is needed to determine the face orientation
const labelList& owner = mesh.faceOwner();
auto cellTypes =
nonNullSmartPointer<vtkUnsignedCharArray>(vtkmesh->GetCellTypesArray());
auto cells =
nonNullSmartPointer<vtkCellArray>(vtkmesh->GetCells());
auto faces =
nonNullSmartPointer<vtkIdTypeArray>(vtkmesh->GetFaces());
auto cellLocations =
nonNullSmartPointer<vtkIdTypeArray>(vtkmesh->GetCellLocationsArray());
auto faceLocations =
nonNullSmartPointer<vtkIdTypeArray>(vtkmesh->GetFaceLocations());
UList<uint8_t> cellTypesUL =
vtkUList
(
cellTypes,
sizing.nFieldCells()
);
UList<vtkIdType> cellsUL =
vtkUList
(
cells,
sizing.nFieldCells(),
sizing.sizeInternal(foamVtuSizing::slotType::CELLS)
);
UList<vtkIdType> cellLocationsUL =
vtkUList
(
cellLocations,
sizing.sizeInternal(foamVtuSizing::slotType::CELLS_OFFSETS)
);
UList<vtkIdType> facesUL =
vtkUList
(
faces,
sizing.sizeInternal(foamVtuSizing::slotType::FACES)
);
UList<vtkIdType> faceLocationsUL =
vtkUList
(
faceLocations,
sizing.sizeInternal(foamVtuSizing::slotType::FACES_OFFSETS)
);
sizing.populateInternal
(
mesh,
cellTypesUL,
cellsUL,
cellLocationsUL,
facesUL,
faceLocationsUL,
static_cast<foamVtkMeshMaps&>(vtuData)
);
labelList& cellMap = vtuData.cellMap();
labelList& addPointCellLabels = vtuData.additionalIds();
// Scan for cells which need to be decomposed and count additional points
// and cells
if (!reader_->GetUseVTKPolyhedron())
// Convert OpenFOAM mesh vertices to VTK
// - can only do this *after* populating the decompInfo with cell-ids
// for any additional points (ie, mesh cell-centres)
vtkSmartPointer<vtkPoints> vtkpoints = vtkmesh->GetPoints();
if (!vtkpoints)
{
forAll(cellShapes, celli)
{
const cellModel& model = cellShapes[celli].model();
if
(
model != hex
&& model != wedge
&& model != prism
&& model != pyr
&& model != tet
&& model != tetWedge
)
{
const cell& cFaces = mesh.cells()[celli];
forAll(cFaces, cFacei)
{
const face& f = mesh.faces()[cFaces[cFacei]];
label nQuads = 0;
label nTris = 0;
f.nTrianglesQuads(mesh.points(), nTris, nQuads);
nAddCells += nQuads + nTris;
}
nAddCells--;
nAddPoints++;
}
}
// No points previously, add an entry
vtkpoints = vtkSmartPointer<vtkPoints>::New();
vtkmesh->SetPoints(vtkpoints);
}
// 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> vtkpoints = vtkSmartPointer<vtkPoints>::New();
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<vtkIdType> 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;