Incorrect Volume for QuadraticHexahedron ?!
Hi, it appears that vtk gives out an incorrect volume for QuadraticHexahedrons, i.e. for the Cell defined in vtkQuadraticHexahedron.h.
This is verified by the failure of the following test ‘test_volume.cpp’
#define BOOST_TEST_MODULE test_volume
#include <boost/test/unit_test.hpp>
#include <vtkCellData.h>
#include <vtkIntegrateAttributes.h>
#include <vtkQuadraticHexahedron.h>
#include <vtkQuadraticTetra.h>
#include <vtkUnstructuredGrid.h>
#include <vtkXMLUnstructuredGridWriter.h>
#include <iostream>
const std::string root = PROJECT_ROOT;
double vtk_volume(vtkSmartPointer<vtkUnstructuredGrid> grid)
{
vtkNew<vtkIntegrateAttributes> integrator;
integrator->SetInputData(grid);
integrator->Update();
vtkUnstructuredGrid* output = integrator->GetOutput();
double volume = output->GetCellData()->GetArray("Volume")->GetComponent(0, 0);
return volume;
}
BOOST_AUTO_TEST_CASE(test_quadratic_hex)
{
std::map<int, std::array<double, 3>> nodes = {{486, {80.0, 80.0, 80.0}},
{491, {80.0, 100.0, 80.0}}, {492, {80.0, 90.0, 80.0}},
{508, {100.0, 80.0, 80.0}}, {511, {90.0, 80.0, 80.0}},
{512, {110.0, 100.0, 80.0}}, {513, {100.0, 90.0, 80.0}},
{515, {90.0, 100.0, 80.0}}, {592, {80.0, 80.0, 100.0}},
{596, {80.0, 80.0, 90.0}}, {599, {80.0, 100.0, 100.0}},
{601, {80.0, 90.0, 100.0}}, {603, {80.0, 100.0, 90.0}},
{713, {100.0, 80.0, 100.0}}, {715, {100.0, 80.0, 90.0}},
{716, {90.0, 80.0, 100.0}}, {717, {110.0, 100.0, 100.0}},
{718, {100.0, 90.0, 100.0}}, {719, {110.0, 100.0, 90.0}},
{720, {90.0, 100.0, 100.0}}};
std::vector<int> elem = {
486, 592, 713, 508, 491, 599, 717, 512, // 8 corners
596, 716, 715, 511, 603, 720, 719, // mids
515, 492, 601, 718, 513 // mids
};
auto points = vtkSmartPointer<vtkPoints>::New();
std::map<int, vtkIdType> idMap;
vtkIdType vtkId = 0;
for (auto& p : nodes) {
points->InsertNextPoint(p.second[0], p.second[1], p.second[2]);
idMap[p.first] = vtkId++;
}
auto hex = vtkSmartPointer<vtkQuadraticHexahedron>::New();
for (unsigned int i = 0; i < elem.size(); ++i) {
hex->GetPointIds()->SetId(i, idMap[elem[i]]);
}
auto ug = vtkSmartPointer<vtkUnstructuredGrid>::New();
ug->SetPoints(points);
ug->InsertNextCell(hex->GetCellType(), hex->GetPointIds());
// One calculates by hand that the correct volume is 8666.666
// The method vtk_volume gives out 9333.333
BOOST_CHECK_CLOSE(vtk_volume(ug), 8666.666, 1e-4);
}
By integrating over the quadratic function for the ‘pulled out’ node, one easily verifies that the correct volume is 8666.666!
The same result is obtained by passing this Hexahedron to Calculix (cgx) or Paraview and using their volume functions. Both cgx and paraview give out 8666.666.
Follow this procedure to verify the volume with Calculix (cgx):
-
Consider the following .inp file:
*NODE 486, 80.0 , 80.0 , 80.0 491, 80.0 , 100.0 , 80.0 492, 80.0 , 90.0 , 80.0 508, 100.0 , 80.0 , 80.0 511, 90.0 , 80.0 , 80.0 512, 110.0 , 100.0 , 80.0 513, 100.0 , 90.0 , 80.0 515, 90.0 , 100.0 , 80.0 592, 80.0 , 80.0 , 100.0 596, 80.0 , 80.0 , 90.0 599, 80.0 , 100.0 , 100.0 601, 80.0 , 90.0 , 100.0 603, 80.0 , 100.0 , 90.0 713, 100.0 , 80.0 , 100.0 715, 100.0 , 80.0 , 90.0 716, 90.0 , 80.0 , 100.0 717, 110.0 , 100.0 , 100.0 718, 100.0 , 90.0 , 100.0 719, 110.0 , 100.0 , 90.0 720, 90.0 , 100.0 , 100.0 **HW_COMPONENT ID=1 COLOR=4 NAME=cube *ELEMENT,TYPE=C3D20 140, 486, 592, 713, 508, 491, 599, 717, 512, 596, 716, 715, 511, 603, 720, 719, 515, 492, 601, 718, 513 **HW_SET ID=5 *NSET, NSET=NFIX_Z **HW_SET ID=4 *NSET, NSET=NFIX_Y **HW_SET ID=3 *NSET, NSET=NFIX_X **HW_SET ID=2 *NSET, NSET=NALL 486, 491, 492, 508, 511, 512, 513, 515, 592, 596, 599, 601, 603, 713, 715, 716, 717, 718, 719, 720 *ELSET, ELSET=cube 140, -
You may save this .inp file as ‘hex.inp’.
-
Open the file in cgx:
cgx -c ‘hex.inp’ -
In the cgx-Viewer left-click and hold to select
Toggle-Commandline. You should see:at the bottom of the window now. This is the cgx-commandline -
At the cgx-commandline type
volu cubeand hit enter -
At the window you launched cgx from you should now see the output
Elements: 1 Nodes:20 Datasets:0 MinElemNr: 140 MaxElemNr: 140 MinNodNr:486 MaxNodNr:720 read in 0.020798 sec found elements of type 1: 0 found elements of type 2: 0 found elements of type 3: 0 found elements of type 4: 1 found elements of type 5: 0 found elements of type 6: 0 found elements of type 7: 0 found elements of type 8: 0 found elements of type 9: 0 found elements of type 10: 0 found elements of type 11: 0 found elements of type 12: 0 gtol calculated:3.000000e-04 VOLUME:8.666667e+03 CENTER OF GRAVITY: 9.107692e+01 9.076923e+01 9.000000e+01
The last line indicates that the volume is 8666.666…