Dear all,
This is driving me round the bend. When I have a meshing coming from a
rectilinear source, I get NANs from my normal and face area calculations.
When all other input meshes are used, i.e. unstructured, the results are
correct. To test this I've written the following test code which makes a
single cell mesh and produces the NANs. I then demonstrate that by using
the transform classes, with zero rotation, I can then extract the same mesh
back as a structured grid, and the errors disappear. I can also apply a
non-zero rotation and the errors disappear too.
Note I apply a rotation, zero or otherwise, because my hypothesis having
stepped through 6.1 in the debugger, is that the vtk_pixel cell type (axis
aligned) coming from the rectilinear grid does not play well with the
normal calculation. But I can't find the bug, only that all entries of the
normal are zero and the rest of the function fails.
Clearly, and more likely, is that I could also not be using the classes
correctly, and if so please let me know?? Hope somebody can find the bug,
I've got nowhere on this after a day of looking. Again, anything non-axis
aligned and or non vtk_pixel would appear to work as expected.
Help appreciated,
Cheers again,
Andy
cmake_minimum_required(VERSION 2.8)
PROJECT(weird)
set(CMAKE_COLOR_MAKEFILE ON)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wextra -Wall -Wconversion -pedantic
-Wuninitialized -std=c++11 -Wno-extra-semi")
FIND_PACKAGE(VTK REQUIRED NO_MODULE)
INCLUDE(${VTK_USE_FILE})
find_package(Qt5Widgets REQUIRED)
set(SOURCES main.cxx)
ADD_EXECUTABLE(weird ${SOURCES})
TARGET_LINK_LIBRARIES(weird ${VTK_LIBRARIES})
#include "vtkRectilinearGrid.h"
#include "vtkXMLRectilinearGridWriter.h"
#include "vtkXMLStructuredGridWriter.h"
#include "vtkDoubleArray.h"
#include "vtkSmartPointer.h"
#include "vtkCell.h"
#include "vtkPolygon.h"
#include "vtkPoints.h"
#include "vtkIdTypeArray.h"
#include "vtkTransform.h"
#include "vtkTransformFilter.h"
#include "vtkPointSet.h"
#include "vtkStructuredGrid.h"
//
#include <array>
int main()
{
vtkSmartPointer<vtkRectilinearGrid>rg3D = vtkSmartPointer<vtkRectilinearGrid>::New();
vtkSmartPointer<vtkDoubleArray> xdata = vtkSmartPointer<vtkDoubleArray>::New();
vtkSmartPointer<vtkDoubleArray> ydata = vtkSmartPointer<vtkDoubleArray>::New();
vtkSmartPointer<vtkDoubleArray> zdata = vtkSmartPointer<vtkDoubleArray>::New();
const double deltaX = 0.01;
xdata->InsertNextValue(0.0);
xdata->InsertNextValue(deltaX);
ydata->InsertNextValue(0.0);
ydata->InsertNextValue(deltaX);
zdata->InsertNextValue(0.0);
zdata->InsertNextValue(deltaX);
rg3D->SetDimensions(2,2,2);
rg3D->SetXCoordinates(xdata);
rg3D->SetYCoordinates(ydata);
rg3D->SetZCoordinates(zdata);
vtkSmartPointer<vtkXMLRectilinearGridWriter> writer =
vtkSmartPointer<vtkXMLRectilinearGridWriter>::New();
writer->SetInputData(rg3D);
writer->SetFileName("rgrid.vtr");
writer->Write();
vtkSmartPointer<vtkTransform> transform =
vtkSmartPointer<vtkTransform>::New();
//transform->RotateZ(1.2345); // makes no difference
transform->RotateZ(0.0);
vtkSmartPointer<vtkTransformFilter> transformFilter =
vtkSmartPointer<vtkTransformFilter>::New();
transformFilter->SetInputData(rg3D);
transformFilter->SetTransform(transform);
transformFilter->Update();
vtkSmartPointer<vtkXMLStructuredGridWriter> writer2 =
vtkSmartPointer<vtkXMLStructuredGridWriter>::New();
writer2->SetInputData(transformFilter->GetStructuredGridOutput());
writer2->SetFileName("sgrid.vts");
writer2->Write();
/////// switch these two round to enable weirdness.
//auto grid = transformFilter->GetStructuredGridOutput(); // removes NaNs
auto grid = rg3D; // gives NaNs
////////
for(uint c = 0; c < grid->GetNumberOfCells(); c++)
{
auto cell = grid->GetCell(c);
const auto numFaces = cell->GetNumberOfFaces();
for(uint f = 0; f < numFaces; f++)
{
const auto face = cell->GetFace(f);
//std::cout << cell->GetCellType() << " " << face->GetCellType() << std::endl;
vtkSmartPointer<vtkPolygon> polygon = vtkSmartPointer<vtkPolygon>::New();
polygon->GetPoints()->SetNumberOfPoints(4);
polygon->GetPointIds()->SetNumberOfIds(4);
for(uint i = 0; i < face->GetNumberOfPoints(); i++)
{
polygon->GetPointIds()->SetId(i,i);
}
for(uint p = 0; p < face->GetNumberOfPoints(); p++)
{
std::array<double,3> pp;
face->GetPoints()->GetPoint(p, pp.data());
polygon->GetPoints()->SetPoint(p, pp[0], pp[1], pp[2]);
}
const double area = polygon->ComputeArea();
std::cout << "Area " << area << std::endl;
std::array<double,3> faceNormal;
const double faceArea = vtkPolygon::ComputeArea(polygon->GetPoints(),
face->GetNumberOfPoints(),
polygon->GetPointIds()->GetPointer(0),
faceNormal.data());
std::cout << "Area 2 " << faceArea << std::endl;
std::cout << "Face Normal " << faceNormal[0] << " " << faceNormal[1] << " " << faceNormal[2] << std::endl;
vtkSmartPointer<vtkIdTypeArray> pointsIDs = vtkSmartPointer<vtkIdTypeArray>::New();
for(uint i = 0; i < face->GetNumberOfPoints(); i++)
{
pointsIDs->InsertNextValue(i);
}
std::array<double,3> faceCentre;
vtkPolygon::ComputeCentroid(pointsIDs,
polygon->GetPoints(),
faceCentre.data());
std::cout << "Face Centre " << faceCentre[0] << " " << faceCentre[1] << " " << faceCentre[2] << std::endl;
}
}
}
_______________________________________________
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 ParaView Wiki at:
http://paraview.org/Wiki/ParaView
Follow this link to subscribe/unsubscribe:
http://public.kitware.com/mailman/listinfo/paraview
