Commit 54f52ce1 authored by Joachim Pouderoux's avatar Joachim Pouderoux
Browse files

Enhance and optimize Xdmf2 reader and writer to support static meshes.

Reader now caches previous timestep dataset geometry and topology and do not
read them again if new dataset use the same heavy datapath.

Writer has a new property MeshStaticOverTime to specify that the mesh
is static. If so, topology and geometry heavy data will be produced once at t=0
and will be reused for all time steps.
parent 922d08c0
......@@ -43,6 +43,8 @@
#include <deque>
#include <cassert>
#include <libxml/tree.h>
using namespace xdmf2;
static void vtkScaleExtents(int in_exts[6], int out_exts[6], int stride[3])
......@@ -137,7 +139,7 @@ vtkDataObject* vtkXdmfHeavyData::ReadData()
}
//----------------------------------------------------------------------------
vtkDataObject* vtkXdmfHeavyData::ReadData(XdmfGrid* xmfGrid)
vtkDataObject* vtkXdmfHeavyData::ReadData(XdmfGrid* xmfGrid, int blockId)
{
if (!xmfGrid || xmfGrid->GetGridType() == XDMF_GRID_UNSET)
{
......@@ -151,7 +153,7 @@ vtkDataObject* vtkXdmfHeavyData::ReadData(XdmfGrid* xmfGrid)
{
// grid is a temporal collection, pick the sub-grid with matching time and
// process that.
return this->ReadTemporalCollection(xmfGrid);
return this->ReadTemporalCollection(xmfGrid, blockId);
}
else if (gridType == XDMF_GRID_COLLECTION ||
gridType == XDMF_GRID_TREE)
......@@ -160,7 +162,7 @@ vtkDataObject* vtkXdmfHeavyData::ReadData(XdmfGrid* xmfGrid)
}
// grid is a primitive grid, so read the data.
return this->ReadUniformData(xmfGrid);
return this->ReadUniformData(xmfGrid, blockId);
}
//----------------------------------------------------------------------------
......@@ -189,7 +191,7 @@ vtkDataObject* vtkXdmfHeavyData::ReadComposite(XdmfGrid* xmfComposite)
if (!child_is_leaf || !distribute_leaf_nodes ||
(number_of_leaf_nodes % this->NumberOfPieces) == this->Piece)
{
vtkDataObject* childDO = this->ReadData(xmfChild);
vtkDataObject* childDO = this->ReadData(xmfChild, cc);
if (childDO)
{
multiBlock->SetBlock(cc, childDO);
......@@ -204,7 +206,7 @@ vtkDataObject* vtkXdmfHeavyData::ReadComposite(XdmfGrid* xmfComposite)
//----------------------------------------------------------------------------
vtkDataObject* vtkXdmfHeavyData::ReadTemporalCollection(
XdmfGrid* xmfTemporalCollection)
XdmfGrid* xmfTemporalCollection, int blockId)
{
assert(xmfTemporalCollection->GetGridType() & XDMF_GRID_COLLECTION &&
xmfTemporalCollection->GetCollectionType() == XDMF_GRID_COLLECTION_TEMPORAL
......@@ -254,7 +256,7 @@ vtkDataObject* vtkXdmfHeavyData::ReadTemporalCollection(
std::deque<XdmfGrid*>::iterator iter;
for (iter = valid_children.begin(); iter != valid_children.end(); ++iter)
{
vtkDataObject* childDO = this->ReadData(*iter);
vtkDataObject* childDO = this->ReadData(*iter, blockId);
if (childDO)
{
child_data_objects.push_back(childDO);
......@@ -286,7 +288,7 @@ vtkDataObject* vtkXdmfHeavyData::ReadTemporalCollection(
//----------------------------------------------------------------------------
// Read a non-composite grid. Note here uniform has nothing to do with
// vtkUniformGrid but to what Xdmf's GridType="Uniform".
vtkDataObject* vtkXdmfHeavyData::ReadUniformData(XdmfGrid* xmfGrid)
vtkDataObject* vtkXdmfHeavyData::ReadUniformData(XdmfGrid* xmfGrid, int blockId)
{
assert(xmfGrid->IsUniform() && "Input must be a uniform xdmf grid.");
......@@ -300,6 +302,79 @@ vtkDataObject* vtkXdmfHeavyData::ReadUniformData(XdmfGrid* xmfGrid)
// Read heavy data for grid geometry/topology. This does not read any
// data-arrays. They are read explicitly.
XdmfTopology* topo = xmfGrid->GetTopology();
XdmfGeometry* geom = xmfGrid->GetGeometry();
xmlChar* filePtr;
bool caching = true;
XdmfDOM* topoDom = topo->GetDOM();
XdmfXmlNode topoNode = topo->GetElement();
XdmfXmlNode topoNodeDataItem = topoDom->FindElement("DataItem", 0, topoNode);
std::string topoFilename = "NULL";
if (topoNodeDataItem && caching)
{
filePtr = topoNodeDataItem->children->content;
if (filePtr != NULL)
{
topoFilename = reinterpret_cast<char*>(filePtr);
}
else
{
vtkErrorWithObjectMacro(this->Reader, << "Cannot find DataItem element in topology xml, no caching possible");
caching = false;
}
}
else
{
caching = false;
}
XdmfDOM* geomDom = geom->GetDOM();
XdmfXmlNode geomNode = geom->GetElement();
XdmfXmlNode geomNodeDataItem = geomDom->FindElement("DataItem", 0, geomNode);
std::string geomFilename = "NULL";
if (geomNodeDataItem && caching)
{
filePtr = geomNodeDataItem->children->content;
if (filePtr != NULL)
{
geomFilename = reinterpret_cast<char*>(filePtr);
}
else
{
vtkErrorWithObjectMacro(this->Reader, << "Cannot find DataItem element in geometry xml, no caching possible");
caching = false;
}
}
else
{
caching = false;
}
vtkXdmfReader::XdmfReaderCachedData& cache =
vtkXdmfReader::SafeDownCast(this->Reader)->GetDataSetCache();
vtkXdmfReader::XdmfDataSetTopoGeoPath cachedData = cache[blockId];
if (caching &&
(cachedData.topologyPath == topoFilename) && (cachedData.geometryPath == geomFilename))
{
vtkDataSet* ds = vtkDataSet::SafeDownCast(
vtkDataObjectTypes::NewDataObject(cachedData.dataset->GetDataObjectType()));
ds->ShallowCopy(cachedData.dataset);
this->ReadAttributes(ds, xmfGrid);
return ds;
}
if (caching)
{
cachedData.topologyPath = topoFilename;
cachedData.geometryPath = geomFilename;
if (cache[blockId].dataset != NULL)
{
cache[blockId].dataset->Delete();
cache[blockId].dataset = NULL;
}
}
XdmfInt32 status = xmfGrid->Update();
if (status == XDMF_FAIL)
{
......@@ -310,31 +385,36 @@ vtkDataObject* vtkXdmfHeavyData::ReadUniformData(XdmfGrid* xmfGrid)
switch (vtk_data_type)
{
case VTK_UNIFORM_GRID:
dataObject = this->RequestImageData(xmfGrid, true);
break;
case VTK_UNIFORM_GRID:
dataObject = this->RequestImageData(xmfGrid, true);
break;
case VTK_IMAGE_DATA:
dataObject = this->RequestImageData(xmfGrid, false);
break;
case VTK_IMAGE_DATA:
dataObject = this->RequestImageData(xmfGrid, false);
break;
case VTK_STRUCTURED_GRID:
dataObject = this->RequestStructuredGrid(xmfGrid);
break;
case VTK_STRUCTURED_GRID:
dataObject = this->RequestStructuredGrid(xmfGrid);
break;
case VTK_RECTILINEAR_GRID:
dataObject = this->RequestRectilinearGrid(xmfGrid);
break;
case VTK_RECTILINEAR_GRID:
dataObject = this->RequestRectilinearGrid(xmfGrid);
break;
case VTK_UNSTRUCTURED_GRID:
dataObject = this->ReadUnstructuredGrid(xmfGrid);
break;
case VTK_UNSTRUCTURED_GRID:
dataObject = this->ReadUnstructuredGrid(xmfGrid);
break;
default:
// un-handled case.
return 0;
default:
// un-handled case.
return 0;
}
if (caching)
{
cache[blockId].dataset = vtkDataSet::SafeDownCast(dataObject);
dataObject->Register(0);
}
return dataObject;
}
......
......@@ -55,7 +55,7 @@ public:
~vtkXdmfHeavyData();
// Description:
vtkDataObject* ReadData(xdmf2::XdmfGrid* xmfGrid);
vtkDataObject* ReadData(xdmf2::XdmfGrid* xmfGrid, int blockId = -1);
// Description:
vtkDataObject* ReadData();
......@@ -70,10 +70,11 @@ public:
// of points possible.
static int GetNumberOfPointsPerCell(int vtk_cell_type);
private:
// Description:
// Read a temporal collection.
vtkDataObject* ReadTemporalCollection(xdmf2::XdmfGrid* xmfTemporalCollection);
vtkDataObject* ReadTemporalCollection(xdmf2::XdmfGrid* xmfTemporalCollection, int blockId);
// Description:
// Read a spatial-collection or a tree.
......@@ -82,7 +83,7 @@ private:
// Description:
// Read a non-composite grid. Note here uniform has nothing to do with
// vtkUniformGrid but to what Xdmf's GridType="Uniform".
vtkDataObject* ReadUniformData(xdmf2::XdmfGrid* xmfGrid);
vtkDataObject* ReadUniformData(xdmf2::XdmfGrid* xmfGrid, int blockId);
// Description:
// Reads the topology and geometry for an unstructured grid. Does not read any
......@@ -156,7 +157,6 @@ private:
// Creates a new dataset with egdes selected by the set, extracting them from
// the input dataset.
vtkDataSet* ExtractEdges(xdmf2::XdmfSet* xmfSet, vtkDataSet* dataSet);
};
#endif
......@@ -19,6 +19,7 @@
#include "vtkCharArray.h"
#include "vtkCompositeDataPipeline.h"
#include "vtkDataObjectTypes.h"
#include "vtkDataSet.h"
#include "vtkExtentTranslator.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
......@@ -116,6 +117,8 @@ vtkXdmfReader::~vtkXdmfReader()
delete this->CellArraysCache;
delete this->GridsCache;
delete this->SetsCache;
this->ClearDataSetCache();
}
//----------------------------------------------------------------------------
......@@ -318,9 +321,15 @@ int vtkXdmfReader::RequestInformation(vtkInformation *, vtkInformationVector **,
// * Publish the SIL which provides information about the grid hierarchy.
outInfo->Set(vtkDataObject::SIL(), domain->GetSIL());
const std::map<int, XdmfFloat64>& ts = domain->GetTimeStepsRev();
std::vector<double> time_steps(ts.size());
std::map<int, XdmfFloat64>::const_iterator it = ts.begin();
for (int i = 0; it != ts.end(); i++, ++it)
time_steps[i] = it->second;
// * Publish time information.
std::vector<double> time_steps(domain->GetTimeSteps().begin(),
domain->GetTimeSteps().end());
/*std::vector<double> time_steps(domain->GetTimeSteps().begin(),
domain->GetTimeSteps().end());*/
if (time_steps.size() > 0)
{
......@@ -386,6 +395,10 @@ int vtkXdmfReader::RequestData(vtkInformation *, vtkInformationVector **,
}
this->LastTimeIndex = this->ChooseTimeStep(outInfo);
if (this->LastTimeIndex == 0)
{
this->ClearDataSetCache();
}
vtkXdmfHeavyData dataReader(this->XdmfDocument->GetActiveDomain(), this);
dataReader.Piece = updatePiece;
......@@ -615,3 +628,25 @@ vtkGraph* vtkXdmfReader::GetSIL()
}
return 0;
}
//----------------------------------------------------------------------------
void vtkXdmfReader::ClearDataSetCache()
{
XdmfReaderCachedData::iterator it = this->DataSetCache.begin();
while (it != this->DataSetCache.end())
{
if (it->second.dataset != NULL)
{
it->second.dataset->Delete();
}
++it;
}
this->DataSetCache.clear();
}
//----------------------------------------------------------------------------
vtkXdmfReader::XdmfReaderCachedData& vtkXdmfReader::GetDataSetCache()
{
return this->DataSetCache;
}
......@@ -35,6 +35,7 @@
#include "vtkIOXdmf2Module.h" // For export macro
#include "vtkDataReader.h"
#include <map> // for caching
class vtkXdmfArraySelection;
class vtkXdmfDocument;
......@@ -140,6 +141,21 @@ public:
virtual vtkGraph* GetSIL();
//BTX
class XdmfDataSetTopoGeoPath
{
public:
XdmfDataSetTopoGeoPath() : dataset(0), topologyPath(), geometryPath() {}
vtkDataSet* dataset;
std::string topologyPath;
std::string geometryPath;
};
typedef std::map<int, XdmfDataSetTopoGeoPath> XdmfReaderCachedData;
// Description
// Get the data set cache
XdmfReaderCachedData& GetDataSetCache();
protected:
vtkXdmfReader();
~vtkXdmfReader();
......@@ -182,11 +198,16 @@ protected:
vtkXdmfArraySelection* SetsCache;
int SILUpdateStamp;
XdmfReaderCachedData DataSetCache;
private:
// Description:
// Prepares the XdmfDocument.
bool PrepareDocument();
void ClearDataSetCache();
// Description:
// Returns the time-step index requested using the UPDATE_TIME_STEPS from the
// information.
......
......@@ -325,7 +325,13 @@ int vtkXdmfDomain::GetVTKDataType(XdmfGrid* xmfGrid)
//----------------------------------------------------------------------------
int vtkXdmfDomain::GetIndexForTime(double time)
{
std::set<XdmfFloat64>::iterator iter = this->TimeSteps.upper_bound(time);
std::map<XdmfFloat64, int>::const_iterator iter = this->TimeSteps.find(time);
if (iter != this->TimeSteps.end())
{
return iter->second;
}
iter = this->TimeSteps.upper_bound(time);
if (iter == this->TimeSteps.begin())
{
// The requested time step is before any available time. We will use it by
......@@ -337,7 +343,7 @@ int vtkXdmfDomain::GetIndexForTime(double time)
iter--;
}
std::set<XdmfFloat64>::iterator iter2 = this->TimeSteps.begin();
std::map<XdmfFloat64, int>::iterator iter2 = this->TimeSteps.begin();
int counter = 0;
while (iter2 != iter)
{
......@@ -641,7 +647,12 @@ void vtkXdmfDomain::CollectNonLeafMetaData(XdmfGrid* xmfGrid,
XdmfTime* xmfTime = xmfGrid->GetTime();
if (xmfTime && xmfTime->GetTimeType() != XDMF_TIME_UNSET)
{
this->TimeSteps.insert(xmfTime->GetValue());
int step = static_cast<int>(this->TimeSteps.size());
if (this->TimeSteps.find(xmfTime->GetValue()) == this->TimeSteps.end())
{
this->TimeSteps[xmfTime->GetValue()] = step;//this->TimeSteps.insert(xmfTime->GetValue());
this->TimeStepsRev[step] = xmfTime->GetValue();
}
}
}
}
......@@ -728,7 +739,12 @@ void vtkXdmfDomain::CollectLeafMetaData(XdmfGrid* xmfGrid, vtkIdType silParent)
XdmfTime* xmfTime = xmfGrid->GetTime();
if (xmfTime && xmfTime->GetTimeType() != XDMF_TIME_UNSET)
{
this->TimeSteps.insert(xmfTime->GetValue());
int step = static_cast<int>(this->TimeSteps.size());
if (this->TimeSteps.find(xmfTime->GetValue()) == this->TimeSteps.end())
{
this->TimeSteps[xmfTime->GetValue()] = step;//this->TimeSteps.insert(xmfTime->GetValue());
this->TimeStepsRev[step] = xmfTime->GetValue();
}
}
}
......
......@@ -203,8 +203,9 @@ private:
vtkXdmfArraySelection* CellArrays;
vtkXdmfArraySelection* Grids;
vtkXdmfArraySelection* Sets;
std::set<XdmfFloat64> TimeSteps; //< Only discrete timesteps are currently
std::map<XdmfFloat64, int> TimeSteps; //< Only discrete timesteps are currently
// supported.
std::map<int, XdmfFloat64> TimeStepsRev;
public:
//---------------------------------------------------------------------------
......@@ -244,8 +245,10 @@ public:
//---------------------------------------------------------------------------
// Description:
// Returns the timesteps.
const std::set<XdmfFloat64>& GetTimeSteps()
const std::map<XdmfFloat64, int>& GetTimeSteps()
{ return this->TimeSteps; }
const std::map<int,XdmfFloat64>& GetTimeStepsRev()
{ return this->TimeStepsRev; }
//---------------------------------------------------------------------------
// Description:
......@@ -257,8 +260,8 @@ public:
// Returns the time value at the given index.
XdmfFloat64 GetTimeForIndex(int index)
{
std::set<XdmfFloat64>::iterator iter;
int cc=0;
std::map<int, XdmfFloat64>::iterator iter = this->TimeStepsRev.find(index);
/*int cc=0;
for (iter = this->TimeSteps.begin(); iter != this->TimeSteps.end();
iter++, cc++)
{
......@@ -266,7 +269,9 @@ public:
{
return *iter;
}
}
}*/
if (iter != this->TimeStepsRev.end())
return iter->second;
// invalid index.
return 0.0;
}
......
......@@ -53,9 +53,12 @@
#include "XdmfTime.h"
#include "XdmfTopology.h"
#include "vtksys/SystemTools.hxx"
#include <algorithm>
#include <map>
#include <stdio.h>
#include <sstream>
#include <vector>
#include <libxml/tree.h> // always after std::blah stuff
......@@ -186,6 +189,7 @@ vtkXdmfWriter::vtkXdmfWriter()
this->TopTemporalGrid = NULL;
this->DomainMemoryHandler = NULL;
this->SetNumberOfOutputPorts(0);
this->MeshStaticOverTime = false;
}
//----------------------------------------------------------------------------
......@@ -248,6 +252,10 @@ int vtkXdmfWriter::Write()
// always write even if the data hasn't changed
this->Modified();
this->TopologyAtT0.clear();
this->GeometryAtT0.clear();
this->UnlabelledDataArrayId = 0;
//TODO: Specify name of heavy data companion file?
if (!this->DOM)
{
......@@ -331,6 +339,15 @@ int vtkXdmfWriter::RequestData(
return 1;
}
this->WorkingDirectory = vtksys::SystemTools::GetFilenamePath(this->FileName);
this->BaseFileName = vtksys::SystemTools::GetFilenameWithoutLastExtension(this->FileName);
// If mesh is static we force heavy data to be exported in HDF
int lightDataLimit = this->LightDataLimit;
this->LightDataLimit = this->MeshStaticOverTime ? 1 : this->LightDataLimit;
this->CurrentBlockIndex = 0;
if (this->CurrentTimeIndex == 0 &&
this->WriteAllTimeSteps &&
this->NumberOfTimeSteps > 1)
......@@ -346,13 +363,13 @@ int vtkXdmfWriter::RequestData(
tgrid->SetDeleteOnGridDelete(true);
tgrid->SetGridType(XDMF_GRID_COLLECTION);
tgrid->SetCollectionType(XDMF_GRID_COLLECTION_TEMPORAL);
tgrid->SetName(this->BaseFileName.c_str());
XdmfTopology *t = tgrid->GetTopology();
t->SetTopologyType(XDMF_NOTOPOLOGY);
XdmfGeometry *geo = tgrid->GetGeometry();
geo->SetGeometryType(XDMF_GEOMETRY_NONE);
this->DomainMemoryHandler->InsertGrid(tgrid);
this->TopTemporalGrid = tgrid;
}
......@@ -367,19 +384,21 @@ int vtkXdmfWriter::RequestData(
this->DomainMemoryHandler->InsertGrid(grid);
}
this->CurrentTime = 0;
vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
vtkDataObject* input = inInfo->Get(vtkDataObject::DATA_OBJECT());
vtkInformation *inDataInfo = input->GetInformation();
if (inDataInfo->Has(vtkDataObject::DATA_TIME_STEP()))
{
//I am assuming we are not given a temporal data object and getting just one time.
double dataT = input->GetInformation()->Get(vtkDataObject::DATA_TIME_STEP());
//cerr << "Writing " << this->CurrentTimeIndex << " " << *dataT << endl;
this->CurrentTime = input->GetInformation()->Get(vtkDataObject::DATA_TIME_STEP());
//cerr << "Writing timestep" << this->CurrentTimeIndex << " (" << this->CurrentTime << ")" << endl;
XdmfTime *xT = grid->GetTime();
xT->SetDeleteOnGridDelete(true);
xT->SetTimeType(XDMF_TIME_SINGLE);
xT->SetValue(dataT);
xT->SetValue(this->CurrentTime);
grid->Insert(xT);
}
......@@ -397,6 +416,8 @@ int vtkXdmfWriter::RequestData(
this->TopTemporalGrid = NULL;
}
this->LightDataLimit = lightDataLimit;
return 1;
}
......@@ -482,6 +503,19 @@ int vtkXdmfWriter::WriteCompositeDataSet(vtkCompositeDataSet *dobj, xdmf2::XdmfG
return 1;
}
//----------------------------------------------------------------------------
void vtkXdmfWriter::SetupDataArrayXML(XdmfElement* e, XdmfArray* a) const
{
std::stringstream ss;
ss << "<DataItem Dimensions = \"" << a->GetShapeAsString() <<
"\" NumberType = \"" << XdmfTypeToClassString(a->GetNumberType()) <<
"\" Precision = \"" << a->GetElementSize() <<
"\" Format = \"HDF\">" <<
a->GetHeavyDataSetName() << "</DataItem>";
e->SetDataXml(ss.str().c_str());
}
//------------------------------------------------------------------------------
int vtkXdmfWriter::CreateTopology(vtkDataSet *ds, xdmf2::XdmfGrid *grid, vtkIdType PDims[3], vtkIdType CDims[3], vtkIdType &PRank, vtkIdType &CRank, void *staticdata)
{
......@@ -494,14 +528,35 @@ int vtkXdmfWriter::CreateTopology(vtkDataSet *ds, xdmf2::XdmfGrid *grid, vtkIdTy
if (this->HeavyDataFileName)
{
heavyDataSetName = std::string(this->HeavyDataFileName) + ":";
if (this->HeavyDataGroupName)
if (this->MeshStaticOverTime)
{
std::stringstream hdf5group;
hdf5group << "/Topology_";
if (this->CurrentBlockIndex >= 0)
{
if (grid->GetName())
{
hdf5group << grid->GetName();
}
else
{
hdf5group << "Block_" << this->CurrentBlockIndex;
}
heavyDataSetName = heavyDataSetName + hdf5group.str();
}
}
else
{
heavyDataSetName = heavyDataSetName + HeavyDataGroupName + "/Topology";
if (this->HeavyDataGroupName)
{
heavyDataSetName = heavyDataSetName + HeavyDataGroupName + "/Topology";
}
}
heavyName = heavyDataSetName.c_str();
}
XdmfTopology *t = grid->GetTopology();
t->SetLightDataLimit(this->LightDataLimit);
//
// If the topology is unchanged from last grid written, we can reuse the XML
......@@ -528,6 +583,28 @@ int vtkXdmfWriter::CreateTopology(vtkDataSet *ds, xdmf2::XdmfGrid *grid, vtkIdTy
}
}
if (this->MeshStaticOverTime)
{
if (this->CurrentTimeIndex == 0)
{
// Save current topology node at t0 for next time steps
this->TopologyAtT0.push_back(t);
}
else if (static_cast<int>(this->TopologyAtT0.size()) > this->CurrentBlockIndex)
{
// Get topology node at t0
XdmfTopology* topo = this->TopologyAtT0[this->CurrentBlockIndex];
// Setup current topology node with t0 properties
t->SetTopologyTypeFromString(topo->GetTopologyTypeAsString());
t->SetNumberOfElements(topo->GetNumberOfElements());
// Setup connectivity data XML according t0 one