Commit 2b788633 authored by David C. Lonie's avatar David C. Lonie
Browse files

Add VASPTessellationReader for NPT_Z_TESSELLATE.out files.

parent 9838dcfa
......@@ -13,6 +13,7 @@ set(Module_SRCS
vtkProteinRibbonFilter.cxx
vtkSimpleBondPerceiver.cxx
vtkVASPAnimationReader.cxx
vtkVASPTessellationReader.cxx
vtkXYZMolReader2.cxx
)
......
set(TestVASPAnimationReader_ARGS DATA{../Data/VASP/NPT_Z_ANIMATE.out})
set(TestVASPTessellationReader_ARGS DATA{../Data/VASP/NPT_Z_TESSELLATE.out})
vtk_add_test_cxx(${vtk-module}CxxTests tests
TestBallAndStick.cxx
TestPDBBallAndStick.cxx
......@@ -19,6 +21,7 @@ vtk_add_test_cxx(${vtk-module}CxxTests tests
TestProteinRibbon.cxx
TestSimpleBondPerceiver.cxx,NO_VALID
TestVASPAnimationReader.cxx
TestVASPTessellationReader.cxx
TestVDWSpheres.cxx
TestCMLMoleculeReader.cxx
)
......
/*=========================================================================
Program: Visualization Toolkit
Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
All rights reserved.
See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notice for more information.
=========================================================================*/
#include "vtkTestUtilities.h"
#include "vtkRegressionTestImage.h"
#include "vtkActor.h"
#include "vtkCamera.h"
#include "vtkExecutive.h"
#include "vtkDataSetSurfaceFilter.h"
#include "vtkInformation.h"
#include "vtkMolecule.h"
#include "vtkMoleculeMapper.h"
#include "vtkNew.h"
#include "vtkPolyData.h"
#include "vtkPolyDataMapper.h"
#include "vtkProperty.h"
#include "vtkRenderWindowInteractor.h"
#include "vtkRenderWindow.h"
#include "vtkRenderer.h"
#include "vtkStreamingDemandDrivenPipeline.h"
#include "vtkVASPTessellationReader.h"
int TestVASPTessellationReader(int argc, char *argv[])
{
if (argc < 2)
{
std::cerr << "Missing test file argument." << std::endl;
return EXIT_FAILURE;
}
std::string fname(argv[1]);
vtkNew<vtkVASPTessellationReader> reader;
reader->SetFileName(fname.c_str());
reader->UpdateInformation();
vtkInformation *outInfo = reader->GetExecutive()->GetOutputInformation(0);
double *times = outInfo->Get(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
int nTimes = outInfo->Length(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
if (nTimes < 8)
{
std::cerr << "Need at least 8 timesteps, only " << nTimes << " found.\n";
return EXIT_FAILURE;
}
vtkNew<vtkDataSetSurfaceFilter> geomFilter;
geomFilter->SetInputConnection(reader->GetOutputPort(1));
// Show different time steps in each renderer:
vtkNew<vtkRenderer> rens[4];
rens[0]->SetViewport(0.0, 0.5, 0.5, 1.0);
rens[1]->SetViewport(0.5, 0.5, 1.0, 1.0);
rens[2]->SetViewport(0.0, 0.0, 0.5, 0.5);
rens[3]->SetViewport(0.5, 0.0, 1.0, 0.5);
vtkNew<vtkMoleculeMapper> molMappers[4];
vtkNew<vtkActor> molActors[4];
vtkNew<vtkPolyDataMapper> tessMappers[4];
vtkNew<vtkActor> tessActors[4];
vtkNew<vtkRenderWindow> win;
for (size_t i = 0; i < 4; ++i)
{
// Render different timestamps for each:
vtkNew<vtkMolecule> mol;
reader->UpdateTimeStep(times[2 * i]);
mol->ShallowCopy(reader->GetOutput(0));
vtkNew<vtkPolyData> polyData;
geomFilter->UpdateTimeStep(times[2 * i]);
polyData->ShallowCopy(geomFilter->GetOutput(0));
// Rendering setup:
molMappers[i]->SetInputData(mol.Get());
molMappers[i]->UseBallAndStickSettings();
molMappers[i]->RenderLatticeOn();
molActors[i]->SetMapper(molMappers[i].Get());
rens[i]->AddActor(molActors[i].Get());
tessMappers[i]->SetInputData(polyData.Get());
tessMappers[i]->SelectColorArray("Atomic Numbers");
tessMappers[i]->SetLookupTable(molMappers[i]->GetLookupTable());
tessActors[i]->SetMapper(tessMappers[i].Get());
tessActors[i]->GetProperty()->SetOpacity(0.5);
rens[i]->AddActor(tessActors[i].Get());
rens[i]->SetBackground(0.0, 0.0, 0.0);
win->AddRenderer(rens[i].Get());
}
vtkNew<vtkRenderWindowInteractor> iren;
iren->SetRenderWindow(win.GetPointer());
win->SetSize(450,450);
win->Render();
for (size_t i = 0; i < 4; ++i)
{
rens[i]->GetActiveCamera()->Dolly(1.5);
rens[i]->ResetCameraClippingRange();
}
win->Render();
// Finally render the scene and compare the image to a reference image
win->SetMultiSamples(0);
win->GetInteractor()->Initialize();
win->GetInteractor()->Start();
return EXIT_SUCCESS;
}
/*=========================================================================
Program: Visualization Toolkit
Module: vtkVASPTessellationReader.cxx
Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
All rights reserved.
See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notice for more information.
=========================================================================*/
#include "vtkVASPTessellationReader.h"
#include "vtkBoundingBox.h"
#include "vtkCellData.h"
#include "vtkCellType.h"
#include "vtkDataSetAttributes.h"
#include "vtkFloatArray.h"
#include "vtkIdTypeArray.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkMolecule.h"
#include "vtkNew.h"
#include "vtkObjectFactory.h"
#include "vtkPointLocator.h"
#include "vtkPoints.h"
#include "vtkStreamingDemandDrivenPipeline.h"
#include "vtkUnsignedShortArray.h"
#include "vtkUnstructuredGrid.h"
#include "vtkVector.h"
#include "vtkVectorOperators.h"
#include <vtksys/RegularExpression.hxx>
#include <algorithm>
#include <cassert>
#include <cmath>
#include <fstream>
#include <set>
#include <sstream>
#include <vector>
vtkStandardNewMacro(vtkVASPTessellationReader)
typedef vtksys::RegularExpression RegEx;
typedef vtkStreamingDemandDrivenPipeline vtkSDDP;
namespace {
template <typename T>
bool parse(const std::string &str, T &result)
{
if (!str.empty())
{
std::istringstream tmp(str);
tmp >> result;
return !tmp.fail();
}
return false;
}
template <typename T>
bool parseCommaSepList(const std::string &input, std::vector<T> &data)
{
std::istringstream in(input);
for (std::string valStr; std::getline(in, valStr, ',');)
{
T tmp;
if (!parse(valStr, tmp))
{
return false;
}
data.push_back(tmp);
}
return true;
}
// This parses the voronoi points/faces list. The input is expected to be:
// [number of commaSepLists], ([commaSepList]) ([commaSepList]) ...
template <typename T>
bool parseVariableLists(const std::string &input,
std::vector<std::vector<T> > &data,
RegEx *parenExtract)
{
std::istringstream in(input);
size_t nLists;
in >> nLists;
if (in.fail())
{
return false;
}
data.resize(nLists);
std::string parseBuffer = input;
for (size_t i = 0; i < nLists; ++i)
{
if (!parenExtract->find(parseBuffer))
{
return false;
}
if (!parseCommaSepList(parenExtract->match(1), data[i]))
{
return false;
}
// Chop the current match off of the buffer to prepare for the next iter.
parseBuffer = parseBuffer.substr(parenExtract->end());
}
return true;
}
} // end anon namespace
//------------------------------------------------------------------------------
void vtkVASPTessellationReader::PrintSelf(std::ostream &os, vtkIndent indent)
{
this->Superclass::PrintSelf(os, indent);
}
//------------------------------------------------------------------------------
vtkVASPTessellationReader::vtkVASPTessellationReader()
: FileName(NULL),
TimeParser(new RegEx("^ *time *= *([0-9EeDd.+-]+) *$")), // time = (timeVal)
LatticeParser(new RegEx("^ *Rxx *= *([0-9EeDd.+-]+) *," // Rxx
" *Rxy *= *([0-9EeDd.+-]+) *," // Rxy
" *Rxz *= *([0-9EeDd.+-]+) *," // Rxz
" *Ryy *= *([0-9EeDd.+-]+) *," // Ryy
" *Ryz *= *([0-9EeDd.+-]+) *," // Ryz
" *Rzz *= *([0-9EeDd.+-]+) *$" // Rzz
)),
AtomCountParser(new RegEx("^ *Natoms *= *([0-9]+) *$")), // Natoms = (int)),
AtomParser(new RegEx("^ *([0-9]+) *," // Atom index
" *\\(" // Open paren
" *([0-9EeDd.+-]+) *," // X coord
" *([0-9EeDd.+-]+) *," // Y coord
" *([0-9EeDd.+-]+)" // Z coord
" *\\) *," // Close paren
" *([0-9EeDd.+-]+) *$" // Radius
)),
ParenExtract(new RegEx("\\(([^(]+)\\)")) // Extract contents of (...)
{
this->SetNumberOfInputPorts(0);
this->SetNumberOfOutputPorts(2);
}
//------------------------------------------------------------------------------
vtkVASPTessellationReader::~vtkVASPTessellationReader()
{
this->SetFileName(NULL);
delete this->TimeParser;
delete this->LatticeParser;
delete this->AtomCountParser;
delete this->AtomParser;
delete this->ParenExtract;
}
//------------------------------------------------------------------------------
int vtkVASPTessellationReader::RequestData(vtkInformation *,
vtkInformationVector **,
vtkInformationVector *outInfos)
{
vtkInformation *outInfo0 = outInfos->GetInformationObject(0);
vtkInformation *outInfo1 = outInfos->GetInformationObject(1);
vtkMolecule *molecule = vtkMolecule::SafeDownCast(
outInfo0->Get(vtkDataObject::DATA_OBJECT()));
assert(molecule);
vtkUnstructuredGrid *voronoi = vtkUnstructuredGrid::SafeDownCast(
outInfo1->Get(vtkDataObject::DATA_OBJECT()));
assert(voronoi);
std::ifstream in(this->FileName);
if (!in)
{
vtkErrorMacro("Could not open file for reading: "
<< (this->FileName ? this->FileName : ""));
return 1;
}
// Advance to the selected timestep:
size_t stepIdx = this->SelectTimeStepIndex(outInfo0);
double time = 0.;
for (size_t i = 0; i <= stepIdx; ++i) // <= to read the "time=" line
{
if (!this->NextTimeStep(in, time))
{
vtkErrorMacro("Error -- attempting to read timestep #" << (stepIdx + 1)
<< " but encountered a parsing error at timestep #"
<< (i + 1) << ".");
return 1;
}
}
if (this->ReadTimeStep(in, molecule, voronoi))
{
molecule->GetInformation()->Set(vtkDataObject::DATA_TIME_STEP(), time);
voronoi->GetInformation()->Set(vtkDataObject::DATA_TIME_STEP(), time);
}
else
{
molecule->Initialize();
voronoi->Initialize();
}
return 1;
}
//------------------------------------------------------------------------------
int vtkVASPTessellationReader::RequestInformation(
vtkInformation *, vtkInformationVector **,
vtkInformationVector *outInfos)
{
std::ifstream in(this->FileName);
if (!in)
{
vtkErrorMacro("Could not open file for reading: "
<< (this->FileName ? this->FileName : ""));
return 1;
}
// Scan the file for timesteps:
double time;
std::vector<double> times;
double timeRange[2] = { VTK_DOUBLE_MAX, VTK_DOUBLE_MIN };
while (this->NextTimeStep(in, time))
{
times.push_back(time);
timeRange[0] = std::min(timeRange[0], time);
timeRange[1] = std::max(timeRange[1], time);
}
if (!times.empty())
{
for (int port = 0; port < 2; ++port)
{
vtkInformation *outInfo = outInfos->GetInformationObject(port);
outInfo->Set(vtkSDDP::TIME_RANGE(), timeRange, 2);
outInfo->Set(vtkSDDP::TIME_STEPS(), &times[0],
static_cast<int>(times.size()));
}
}
return 1;
}
//------------------------------------------------------------------------------
int vtkVASPTessellationReader::FillOutputPortInformation(int port,
vtkInformation *info)
{
switch (port)
{
case 0:
info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkMolecule");
break;
case 1:
info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkUnstructuredGrid");
break;
default:
return 0;
}
return 1;
}
//------------------------------------------------------------------------------
bool vtkVASPTessellationReader::NextTimeStep(std::istream &in, double &time)
{
std::string line;
while (std::getline(in, line))
{
if (this->TimeParser->find(line))
{
// Parse timestamp:
if (!parse(this->TimeParser->match(1), time))
{
vtkErrorMacro("Error parsing time information from line: " << line);
return false;
}
return true;
}
}
return false;
}
//------------------------------------------------------------------------------
size_t vtkVASPTessellationReader::SelectTimeStepIndex(vtkInformation *info)
{
if (!info->Has(vtkSDDP::TIME_STEPS()) ||
!info->Has(vtkSDDP::UPDATE_TIME_STEP()))
{
return 0;
}
double *times = info->Get(vtkSDDP::TIME_STEPS());
int nTimes = info->Length(vtkSDDP::TIME_STEPS());
double t = info->Get(vtkSDDP::UPDATE_TIME_STEP());
double resultDiff = VTK_DOUBLE_MAX;
size_t result = 0;
for (int i = 0; i < nTimes; ++i)
{
double diff = std::fabs(times[i] - t);
if (diff < resultDiff)
{
resultDiff = diff;
result = static_cast<size_t>(i);
}
}
return result;
}
//------------------------------------------------------------------------------
bool vtkVASPTessellationReader::ReadTimeStep(std::istream &in,
vtkMolecule *molecule,
vtkUnstructuredGrid *voronoi)
{
// Assume the 'time = ...' line has already been read.
std::string line;
// Read the lattice info:
if (!std::getline(in, line))
{
vtkErrorMacro("Unexpected EOF while reading lattice info.");
return false;
}
if (!this->LatticeParser->find(line))
{
vtkErrorMacro("Error parsing lattice info from line: " << line);
return false;
}
vtkVector3d latA(0.);
vtkVector3d latB(0.);
vtkVector3d latC(0.);
vtkVector3d latO(0.);
// Rxx
if (!parse(this->LatticeParser->match(1), latA[0]))
{
vtkErrorMacro("Error parsing Rxx component '"
<< this->LatticeParser->match(1) << "' from line: " << line);
return false;
}
// Rxy
if (!parse(this->LatticeParser->match(2), latB[0]))
{
vtkErrorMacro("Error parsing Rxy component '"
<< this->LatticeParser->match(2) << "' from line: " << line);
return false;
}
// Rxz
if (!parse(this->LatticeParser->match(3), latC[0]))
{
vtkErrorMacro("Error parsing Rxz component '"
<< this->LatticeParser->match(3) << "' from line: " << line);
return false;
}
// Ryy
if (!parse(this->LatticeParser->match(4), latB[1]))
{
vtkErrorMacro("Error parsing Ryy component '"
<< this->LatticeParser->match(4) << "' from line: " << line);
return false;
}
// Ryz
if (!parse(this->LatticeParser->match(5), latC[1]))
{
vtkErrorMacro("Error parsing Ryz component '"
<< this->LatticeParser->match(5) << "' from line: " << line);
return false;
}
// Rzz
if (!parse(this->LatticeParser->match(6), latC[2]))
{
vtkErrorMacro("Error parsing Rzz component '"
<< this->LatticeParser->match(6) << "' from line: " << line);
return false;
}
molecule->SetLattice(latA, latB, latC);
molecule->SetLatticeOrigin(latO);
// Number of atoms:
if (!std::getline(in, line))
{
vtkErrorMacro("Unexpected EOF while parsing number of atoms.");
return false;
}
if (!this->AtomCountParser->find(line))
{
vtkErrorMacro("Error parsing number of atoms from line: " << line);
return false;
}
vtkIdType nAtoms;
if (!parse(this->AtomCountParser->match(1), nAtoms))
{
vtkErrorMacro("Error parsing number atoms '"
<< this->AtomCountParser->match(1) << "'' from line: "
<< line);
return false;
}
// Skip 'Atomic_Numbers =' line:
if (!std::getline(in, line))
{
vtkErrorMacro("Unexpected EOF.");
return false;
}
// Atomic numbers:
std::vector<unsigned short> atomicNumbers;
atomicNumbers.reserve(static_cast<size_t>(nAtoms));
if (!std::getline(in, line))
{
vtkErrorMacro("Unexpected EOF while reading atomic number list.");
return false;
}
if (!parseCommaSepList(line, atomicNumbers))
{
vtkErrorMacro("Error while parsing atomic number list: " << line);
return false;
}
// Initialize the molecule with atoms, setting just the atomic number.
// We'll add positions as we parse them later.
if (static_cast<size_t>(nAtoms) != atomicNumbers.size())
{
vtkErrorMacro("Error: expected " << nAtoms << " atomic numbers, but only "
"parsed " << atomicNumbers.size());
return false;
}
vtkVector3f pos(0.f);
for (size_t i = 0; i < atomicNumbers.size(); ++i)
{
molecule->AppendAtom(atomicNumbers[i], pos);
}
// Now for the actual list of atoms/tessellations
vtkNew<vtkFloatArray> radii;