[prev in list] [next in list] [prev in thread] [next in thread]
List: vtkusers
Subject: [vtkusers] Paraview bug found in pvtu files containing polyhedrons
From: Andrew Parker via vtkusers <vtkusers () vtk ! org>
Date: 2016-02-29 16:03:46
Message-ID: 45EAF588-0B4F-45B9-B4AB-F1EC8EC656FF () googlemail ! com
[Download RAW message or body]
Dear All,
I have recently been having problems reading pvtu files in paraview when those files \
contain arbitrary polyhedrons: this is from the output of a large scale code that \
runs in parallel via mpi. To show the problem, I have written a very small sample \
code, attached, that runs in serial but reproduces the salient issues experienced by \
my larger code when run in parallel. It is clear, via the toggling a boolean, that \
the problem happens when vtk_polyhedrons types are used. There is a short \
description at the top of the bug.cxx file, followed by instructions. I use stock \
versions of paraview here, although self compiled ones show no difference.
In short compile the code, run it, and open the pvtu file in paraview. You should \
see 4 cells. Switch the "bug" bool to true, recompile, re-run it, and reopen the pvtu \
file. The screen will be blank.
The only difference is that cells are added to the unstructured grid explicitly as \
vtk_polyhedrons using a face-stream, rather than as hexs. I cannot do this in the \
real code! The rest of the code is just there to produce the "serial" mesh, and \
threshold this mesh to produce two "parallel-partition" meshes. I use the new ghost \
type framework as this is consistent with the real code, and this runs and compiles \
on OS X and gcc. Likewise it has been build against vtk 6.3 and the bug manifests \
itself in paraview 4.4 and 5.0. Visit 2.10 does not have this bug and can correctly \
open the pvtu containing polyhedrons: this is my current workaround.
Can any of the developers shed light on this? Is a fix known, if so when will it be \
released? Do others have this problem [1, 2, 3] outside of those long reported? Does \
anyone else have workarounds?
Any help really appreciated,
Andy
[1] http://public.kitware.com/pipermail/vtkusers/2015-May/090835.html
[2] https://cmake.org/pipermail/paraview/2012-October/026456.html
[3] https://cmake.org/pipermail/paraview/2015-January/032950.html
["bug.cxx" (bug.cxx)]
//
#include "vtkRectilinearGrid.h"
#include "vtkDoubleArray.h"
#include "vtkStructuredGrid.h"
#include "vtkXMLStructuredGridWriter.h"
#include "vtkSmartPointer.h"
#include "vtkUnstructuredGrid.h"
#include "vtkCell.h"
#include "vtkPolygon.h"
#include "vtkIdTypeArray.h"
#include "vtkCellCenters.h"
#include "vtkIntArray.h"
#include "vtkPolyData.h"
#include "vtkMeshQuality.h"
#include "vtkCellArray.h"
#include "vtkPoints.h"
#include "vtkThreshold.h"
#include "vtkXMLUnstructuredGridWriter.h"
#include "vtkCellData.h"
#include "vtkXMLPUnstructuredGridWriter.h"
#include <cmath>
#include <unordered_map>
//
/*
A 4x1x1 cell unstructured grid is created and split into 2 pieces to simulate a bug \
found in a large parallel code when writing out flow solutions plus halo data via \
mpi. Here we do not use mpi, but mimic the effect of only having part of a grid and \
the halos on each processor as shown below. Each grid has one ghost cell from the \
adjacent grid. Each grid is written to a .vtu file and a .pvtu file is written \
pointing to the two files. The whole mesh is also written out in a .vtu file as well.
Serial mesh view: [ 0 ][ 1 ][ 2 ][ 3 ]
Partition 0: [ 0 ][ 1 ][halo=2]
Partition 1: [halo=1][ 2 ][ 3 ]
When the cells are inserted into the respective vtkUnstructuredGrid with the cell \
type specified as VTK_POLYHEDRON, paraview is unable to show the grid when the .pvtu \
file is loaded and the screen is blank but can show the individual meshes if the .vtu \
files are loaded separately. This is the same problem the large scale code has where \
the bug was uncovered.
However, we found that if the cells were inserted using the cell type they think they \
are (here pure hexs) (e.g. cell->GetCellType()), paraview can show the mesh from the \
.pvtu file and all looks correct. In the large scale code all the cells are \
arbitrary polyhedrons so switching cell-types like this is not an option, but here we \
see VTK_POLYHEDRON does indeed seem to be the cause of the problem.
ViSiT (2.10.0) does not reproduce this error and can be successfully used to load the \
.pvtu file and therefore is our current workaround.
Due to this we surmise the bug is in paraview not vtk. Finally, we note that we are \
using the new ghost type framework recently added and wondered if this is the root of \
the problem. Paraview 4.4 and 5.0 have the same bug.
---HOWTO:
This can all be tested by toggling the below "bug" bool. True gives an error, false \
no error. Re-compiling, and re-running after toggling the bool shows the problem. \
Open the resulting pvtu files in paraview each time. All files are overwritten.
*/
const bool bug = false;
// this function produces the mesh, see main below.
void process(vtkUnstructuredGrid* grid, const std::size_t& ghostID, const \
std::size_t& rank) {
auto* gridIndex = \
vtkIntArray::SafeDownCast(grid->GetCellData()->GetArray("Index"));
vtkSmartPointer<vtkUnstructuredGrid> partitionMesh;
partitionMesh = vtkSmartPointer<vtkUnstructuredGrid>::New();
vtkSmartPointer<vtkPoints> outPoints;
outPoints = vtkSmartPointer<vtkPoints>::New();
partitionMesh->SetPoints(outPoints);
vtkSmartPointer<vtkIntArray> Index = vtkSmartPointer<vtkIntArray>::New();
Index->SetNumberOfComponents(1);
Index->SetName("Index");
partitionMesh->GetCellData()->AddArray(Index);
std::unordered_map<std::size_t, std::size_t> oldToNew;
const auto numCells = grid->GetNumberOfCells();
for(auto c = 0; c < numCells; c++)
{
auto* cell = grid->GetCell(c);
auto* points = cell->GetPoints();
auto* pointIds = cell->GetPointIds();
// obtain array of ids for the writer
std::vector<vtkIdType> theIds;
for(auto p = 0; p < points->GetNumberOfPoints(); p++)
{
const auto currentId = pointIds->GetId(p);
auto iter = oldToNew.find(currentId);
if(iter == oldToNew.end())
{
const auto newPointId = outPoints->InsertNextPoint(points->GetPoint(p));
oldToNew.insert(std::make_pair(currentId, newPointId));
theIds.push_back(newPointId);
}
else
{
theIds.push_back(iter->second);
}
}
// obtain array of faces for the writer
const auto numFaces = cell->GetNumberOfFaces();
vtkSmartPointer<vtkCellArray> faces = vtkSmartPointer<vtkCellArray>::New();
for(auto f = 0; f < numFaces; f++)
{
auto* face = cell->GetFace(f);
std::vector<vtkIdType> facePointIds;
auto ids = face->GetPointIds();
for(auto k = 0; k < ids->GetNumberOfIds(); k++)
{
const auto id = ids->GetId(k);
auto iter = oldToNew.find(id);
if(iter != oldToNew.end())
facePointIds.push_back(iter->second);
else
{
std::exit(EXIT_FAILURE);
}
}
faces->InsertNextCell(facePointIds.size(), facePointIds.data());
}
if(bug == true)
{
partitionMesh->InsertNextCell(VTK_POLYHEDRON, theIds.size(), theIds.data(), \
numFaces, faces->GetPointer()); }
else
{
partitionMesh->InsertNextCell(cell->GetCellType(), theIds.size(), \
theIds.data()); }
Index->InsertNextValue(gridIndex->GetValue(c));
}
// write the ghost cell flags
vtkSmartPointer<vtkUnsignedCharArray> ghostCells;
ghostCells = vtkSmartPointer<vtkUnsignedCharArray>::New();
ghostCells->SetName(vtkDataSetAttributes::GhostArrayName());
ghostCells->SetNumberOfComponents(1);
partitionMesh->GetCellData()->AddArray(ghostCells);
vtkUnsignedCharArray* ghosts = partitionMesh->GetCellGhostArray();
assert(ghosts != nullptr);
auto* index = vtkIntArray::SafeDownCast(partitionMesh->GetCellData()->GetArray("Index"));
for(std::size_t i = 0; i < partitionMesh->GetNumberOfCells(); i++)
{
const auto val = index->GetValue(i);
if(val == ghostID)
{
ghostCells->InsertNextValue(vtkDataSetAttributes::DUPLICATECELL);
}
else
ghostCells->InsertNextValue(0);
}
vtkSmartPointer<vtkXMLPUnstructuredGridWriter> pwriter = \
vtkSmartPointer<vtkXMLPUnstructuredGridWriter>::New(); \
pwriter->SetFileName("part.pvtu"); pwriter->SetNumberOfPieces(2);
pwriter->SetGhostLevel(1);
pwriter->SetStartPiece(rank);
pwriter->SetEndPiece(rank);
pwriter->SetInputData(partitionMesh);
pwriter->SetDataModeToAscii();
// pwriter->SetDataModeToBinary();
// pwriter->SetCompressorTypeToZLib();
pwriter->Write();
}
int main(int, char**)
{
// Make a very simple 3D structured mesh
vtkSmartPointer<vtkStructuredGrid> structuredGrid = \
vtkSmartPointer<vtkStructuredGrid>::New(); vtkSmartPointer<vtkPoints> points = \
vtkSmartPointer<vtkPoints>::New();
const unsigned int numx = 5;
const unsigned int numy = 2;
const unsigned int numz = 2;
const double xstart = -1.0;
const double xend = 1.0;
const double deltaX = std::abs(xend - xstart) / (double(numx - 1));
const double ystart = 0.0;
const double yend = 0.02;
const double deltaY = std::abs(yend - ystart) / (double(numy - 1));
const double zstart = 0.0;
const double zend = 0.02;
const double deltaZ = std::abs(zend - zstart) / (double(numz - 1));
for(unsigned int k = 0; k < numz; k++)
{
const double zpos = zstart + k * deltaZ;
for(unsigned int j = 0; j < numy; j++)
{
const double ypos = ystart + j * deltaY;
for(unsigned int i = 0; i < numx; i++)
{
const double xpos = xstart + i * deltaX;
points->InsertNextPoint(xpos, ypos, zpos);
}
}
}
// Specify the dimensions of the grid
structuredGrid->SetDimensions(numx, numy, numz);
structuredGrid->SetPoints(points);
// Quick cheat to convert to an unstructuredGrid
vtkSmartPointer<vtkIntArray> Index = vtkSmartPointer<vtkIntArray>::New();
Index->SetNumberOfComponents(1);
Index->SetNumberOfTuples(structuredGrid->GetNumberOfCells());
Index->SetName("Index");
for(int i = 0; i < structuredGrid->GetNumberOfCells(); i++)
Index->SetValue(i, i);
structuredGrid->GetCellData()->AddArray(Index);
auto* cellData = structuredGrid->GetCellData();
cellData->SetActiveScalars("Index");
vtkSmartPointer<vtkThreshold> threshold = vtkSmartPointer<vtkThreshold>::New();
threshold->SetInputData(structuredGrid);
threshold->SetInputArrayToProcess(0, 0, 0, vtkDataObject::FIELD_ASSOCIATION_CELLS, \
vtkDataSetAttributes::SCALARS);
threshold->ThresholdByUpper(-1);
threshold->Update();
// At this point we have the "serial" mesh in unstructured grid format.
vtkSmartPointer<vtkUnstructuredGrid> umesh = threshold->GetOutput();
{
// write the serial mesh just for checking
vtkSmartPointer<vtkXMLUnstructuredGridWriter> writer2 = \
vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New(); writer2->SetInputData(umesh);
writer2->SetFileName("serial.vtu");
writer2->SetCompressorTypeToZLib();
writer2->SetDataModeToBinary();
writer2->Write();
}
//
// split the mesh into 2 pieces each with 2 cells + 1 ghost cell
//
{
// first mesh: use threshold to get first 3 cells
vtkSmartPointer<vtkThreshold> threshold = vtkSmartPointer<vtkThreshold>::New();
threshold->SetInputData(umesh);
threshold->SetInputArrayToProcess(0, 0, 0, \
vtkDataObject::FIELD_ASSOCIATION_CELLS, vtkDataSetAttributes::SCALARS); const \
std::size_t ghostID = 2; threshold->ThresholdByLower(ghostID); // note
threshold->Update();
auto* grid = threshold->GetOutput();
process(grid, ghostID, 0);
}
//
// Now do the same but for the second mesh
//
{
// second mesh: use threshold to get last 3 cells
vtkSmartPointer<vtkThreshold> threshold = vtkSmartPointer<vtkThreshold>::New();
threshold->SetInputData(umesh);
threshold->SetInputArrayToProcess(0, 0, 0, \
vtkDataObject::FIELD_ASSOCIATION_CELLS, vtkDataSetAttributes::SCALARS); const \
std::size_t ghostID = 1; threshold->ThresholdByUpper(ghostID); // note
threshold->Update();
auto* grid = threshold->GetOutput();
process(grid, ghostID, 1);
}
}
["CMakeLists.txt" (CMakeLists.txt)]
cmake_minimum_required(VERSION 2.8)
PROJECT(ptvuBug)
set(CMAKE_COLOR_MAKEFILE ON) \
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wextra -Wall -Wconversion -Wuninitialized \
-std=c++11 -Wno-sign-conversion -Werror=ignored-qualifiers \
-Werror=unused-local-typedefs -Werror=return-type -Werror=cast-qual")
find_package(VTK REQUIRED)
include(${VTK_USE_FILE})
ADD_EXECUTABLE(bug bug.cxx)
TARGET_LINK_LIBRARIES(bug ${VTK_LIBRARIES})
_______________________________________________
Powered by www.kitware.com
Visit other Kitware open-source projects at http://www.kitware.com/opensource/opensource.html
Please keep messages on-topic and check the VTK FAQ at: http://www.vtk.org/Wiki/VTK_FAQ
Search the list archives at: http://markmail.org/search/?q=vtkusers
Follow this link to subscribe/unsubscribe:
http://public.kitware.com/mailman/listinfo/vtkusers
[prev in list] [next in list] [prev in thread] [next in thread]
Configure |
About |
News |
Add a list |
Sponsored by KoreLogic