OilfieldToolsNet/IFEM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Added HDF5 to VTF conversion utility + assosiacted fixes in HDF5 output. This now works for linear/stationary problems and nonlinear/time-dependent problems. But not yet for fully-coupled mixed methods.

Browse Source 
git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@1007 e10b68d5-8a6e-419e-a041-bce267b0401d
This commit is contained in:
kmo 2011-05-22 11:25:42 +00:00 committed by Knut Morten Okstad
parent 8459da8f23
commit dcd55dc7de
14 changed files with 492 additions and 219 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
Apps/FiniteDefElasticity/main_NonLinEl.C
View File

@ -282,7 +282,6 @@ int main (int argc, char** argv)
const double epsT = 1.0e-6;
if (dtDump <= 0.0) dtDump = params.stopTime + 1.0;
if (format < 0) dtSave = params.stopTime + 1.0;
double nextDump = params.time.t + dtDump;
double nextSave = params.time.t + dtSave;
@ -302,7 +301,7 @@ int main (int argc, char** argv)
std::cout <<"\nWriting HDF5 file "<< infile <<".hdf5"<< std::endl;
writer = new DataExporter(true);
writer->registerField("u","solution",DataExporter::SIM, skip2nd ? 0 : -1);
writer->registerField("u","displacement",DataExporter::SIM, skip2nd ? 0 : -1);
writer->setFieldValue("u",model,(void*)&simulator.getSolution());
writer->registerWriter(new HDF5Writer(infile));
writer->registerWriter(new XMLWriter(infile));
@ -337,8 +336,9 @@ int main (int argc, char** argv)
if (params.time.t + epsT*params.time.dt > nextSave)
{
// Save solution variables to VTF for visualization
if (!simulator.saveStep(++iStep,params.time.t,n,skip2nd))
return 6;
if (format >= 0)
if (!simulator.saveStep(++iStep,params.time.t,n,skip2nd))
return 6;
// Save solution variables to HDF5
if (writer)

138
Apps/HDF5toVTF/HDF5toVTF.C Normal file
View File

@ -0,0 +1,138 @@
// $Id$
//==============================================================================
//!
//! \file HDF5toVTF.C
//!
//! \date Apr 08 2011
//!
//! \author Arne Morten Kvarving / SINTEF
//!
//! \brief Convert a HDF5 results database to VTF for visualisation.
//!
//==============================================================================
#include "SIM1D.h"
#include "SIM2D.h"
#include "SIM3D.h"
#include "HDF5Writer.h"
#include "XMLWriter.h"
#include "StringUtils.h"
#include <sstream>
#include <stdlib.h>
#include <string.h>
typedef std::map< std::string,std::vector<XMLWriter::Entry> > ProcessList;
int main (int argc, char** argv)
{
int format = 0;
int n[3] = { 5, 5, 5 };
int dims = 3;
char* infile = 0;
char* vtffile = 0;
for (int i = 1; i < argc; i++)
if (!strcmp(argv[i],"-vtf") && i < argc-1)
format = atoi(argv[++i]);
else if (!strcmp(argv[i],"-nviz") && i < argc-1)
n[0] = n[1] = n[2] = atoi(argv[++i]);
else if (!strcmp(argv[i],"-1D"))
dims = 1;
else if (!strcmp(argv[i],"-2D"))
dims = 2;
else if (!infile)
infile = argv[i];
else if (!vtffile)
vtffile = argv[i];
else
std::cerr <<" ** Unknown option ignored: "<< argv[i] << std::endl;
if (!infile) {
std::cout <<"usage: "<< argv[0]
<<" <inputfile> [<vtffile>] [-nviz <nviz>] [-1D|-2D]"<< std::endl;
return 0;
}
else if (!vtffile)
vtffile = infile;
std::cout <<"\n >>> IFEM HDF5 to VTF converter <<<"
<<"\n ==================================\n"
<<"\nInput file: "<< infile
<<"\nOutput file: "<< vtffile
<<"\nNumber of visualization points: "
<< n[0] <<" "<< n[1] << " " << n[2] << std::endl;
HDF5Writer hdf(strtok(infile,"."),true);
XMLWriter xml(infile);
xml.readInfo();
int levels = xml.getLastTimeLevel();
if (levels > 0) SIMinput::msgLevel = 1;
std::cout <<"Reading "<< infile <<": Time levels = "<< levels << std::endl;
const std::vector<XMLWriter::Entry>& entry = xml.getEntries();
std::vector<XMLWriter::Entry>::const_iterator it;
ProcessList processlist;
for (it = entry.begin(); it != entry.end(); ++it)
processlist[it->patchfile].push_back(*it);
ProcessList::const_iterator pit = processlist.begin();
for (int j = 1; pit != processlist.end(); ++pit, ++j)
{
SIMbase* sim;
if (dims == 1)
sim = new SIM1D;
else if (dims == 2)
sim = new SIM2D;
else
sim = new SIM3D;
std::stringstream dummy;
dummy << "PATCHFILE " << pit->first;
sim->parse(const_cast<char*>(dummy.str().c_str()),dummy);
std::stringstream dummy2;
std::string foo(pit->first);
dummy2 << "NODEFILE " << replaceAll(foo,".g2",".gno");
sim->parse(const_cast<char*>(dummy2.str().c_str()),dummy2);
if (!sim->preprocess(std::vector<int>(),false))
return 1;
bool ok = true;
if (processlist.size() > 1) {
std::string temp(strtok(vtffile,"."));
std::stringstream str;
str <<"-"<< j;
temp.append(str.str());
ok = sim->writeGlv(temp.c_str(),n,format);
}
else
ok = sim->writeGlv(vtffile,n,format);
int block = 0;
for (int i = 0; i <= levels && ok; ++i) {
int k = 20;
if (levels > 0) std::cout <<"\nTime level "<< i << std::endl;
for (it = pit->second.begin(); it != pit->second.end() && ok; ++it) {
Vector vec;
std::cout <<"Reading \""<< it->name <<"\""<< std::endl;
ok = hdf.readField(i,it->name,vec,sim,it->components);
if (it->description == "displacement")
ok &= sim->writeGlvS(vec,n,i+1,block,i,1,NULL,10,it->components);
else if (it->components < 2)
ok &= sim->writeGlvS(vec,n,i+1,block,i,2,it->name.c_str(),k++,1);
else
ok &= sim->writeGlvV(vec,it->name.c_str(),n,i+1,block);
}
if (ok)
ok = sim->writeGlvStep(i+1,i);
}
delete sim;
if (!ok) return 2;
}
return 0;
}

11
CMakeLists.txt
View File

@ -129,12 +129,11 @@ IF(NOT WIN32)
ENDIF(NOT WIN32)
# Make the IFEM library
FILE(GLOB_RECURSE IFEM_SRCS ${PROJECT_SOURCE_DIR}/src/*.[Cf]
${PROJECT_SOURCE_DIR}/3rdparty/*.C)
ADD_LIBRARY(IFEM ${IFEM_SRCS})
# Make the Apps
# Make some Apps
FILE(GLOB_RECURSE Poisson_SRCS ${PROJECT_SOURCE_DIR}/Apps/Poisson/*.C)
ADD_EXECUTABLE(Poisson ${Poisson_SRCS})
TARGET_LINK_LIBRARIES(Poisson IFEM ${DEPLIBS})
@ -145,7 +144,13 @@ FILE(GLOB_RECURSE LinEl_SRCS
ADD_EXECUTABLE(LinEl ${LinEl_SRCS})
TARGET_LINK_LIBRARIES(LinEl IFEM ${DEPLIBS})
# Tests
IF(HDF5_LIBRARIES AND VTFWRITER_LIBRARIES)
FILE(GLOB_RECURSE HDF2VTF_SRCS ${PROJECT_SOURCE_DIR}/Apps/HDF5toVTF/*.C)
ADD_EXECUTABLE(HDF5toVTF ${HDF2VTF_SRCS})
TARGET_LINK_LIBRARIES(HDF5toVTF IFEM ${DEPLIBS})
ENDIF(HDF5_LIBRARIES AND VTFWRITER_LIBRARIES)
# Regression tests
FILE(GLOB_RECURSE LINEL_TESTFILES "${PROJECT_SOURCE_DIR}/Apps/LinearElasticity/Test/*.reg")
FOREACH(TESTFILE ${LINEL_TESTFILES})
ADD_TEST(${TESTFILE} ${PROJECT_SOURCE_DIR}/test/regtest.sh "${CMAKE_BINARY_DIR}/${EXECUTABLE_OUTPUT_PATH}/LinEl" "${TESTFILE}")

28
HOWTO
View File

@ -1,17 +1,17 @@
Dette må gjøres for å kompilere denne versjonen:
sudo apt-get install libboost-dev
sudo apt-get install libblas-dev
sudo apt-get install liblapack-dev
sudo apt-get install libarpack2-dev
sudo apt-get install libsuperlu3-dev
sudo apt-get install libboost-dev
sudo apt-get install libpetsc3.0.0-dev
Installer GoToolsCore og GoTrivariate fra GoTools-svn.
Gjør så cmake ., make og sudo make install i disse katalogene
CMake og out-of-tree builds:
CMake støtter debug og release-builds samtidig vi det som kalles
CMake støtter debug og release-builds samtidig via det som kalles
out-of-tree builds. In-tree builds av App'ene forventer at har byggefilene
i en underkatalog med samme navn som bygge-typen. Feks for å bygge Debug
gjør vi:
@ -27,19 +27,21 @@ det ikke fungere.
Flagg av interesse:
Per default lenker vi mot et minimum av bibliotek. Det betyr ingen
PETSc, ingen SuperLU, ingen VTFWriter og ingen SAMG. Disse kan slås på
med opsjoner:
PETSc og ingen SAMG. Disse kan slås på med opsjoner:
-DENABLE_SUPERLU:BOOL=1, -DENABLE_PETSC:BOOL=1, -DENABLE_VTFWRITER:BOOL=1
og -DENABLE_SAMG:BOOL=1.
-DENABLE_PETSC:BOOL=1 og -DENABLE_SAMG:BOOL=1.
Vi bygger kun libIFEM og Apps/Poisson.
Stokes og FiniteDefElasticity har egen CMakeLists.txt. Disse er
satt opp til å bruke in-tree kopi av libIFEM per default, men sjekker
system hvis den ikke finner in-tree. Du kan tvinge system med
-DFORCE_SYSTEM_IFEM:BOOL=1. Merk at in-tree sjekkes både for <root>/<type> og
i <root>.
For å kompilere mot parallel Petsc bruker du -DENABLE_PARALLEL_PETSC:BOOL=1
Ved å spesifisere -DDISABLE_SUPERLU:BOOL=1 kan du slå av SuperLU ligningsløseren.
Ved å spesifisere -DENABLE_SUPERLU_MT:BOOL=1 kan du aktivere multi-threaded SuperLU
istedet for den serielle versjonen.
Ved å spesifisere -DDISABLE_HDF5:BOOL=1 kan du slå av HDF5-støtten.
For å kompilere mot parallel Petsc bruker du -DENABLE_PARALLEL_PETSC:BOOL=1
Vi bygger kun libIFEM, Apps/Poisson og Apps/LinearElasticity.
Stokes og FiniteDefElasticity har egne CMakeLists.txt. Disse
er satt opp til å bruke in-tree kopi av libIFEM per default,
men sjekker system hvis den ikke finner in-tree.
Du kan tvinge system med -DFORCE_SYSTEM_IFEM:BOOL=1.
Merk at in-tree sjekkes både for <root>/<type> og i <root>.

4
src/ASM/ASMbase.C
View File

@ -484,10 +484,10 @@ bool ASMbase::injectNodeVec (const Vector& nodeVec, Vector& globRes,
{
if (nndof == 0) nndof = nf;
if (nodeVec.size() < MLGN.size()*nndof)
if (nodeVec.size() != MLGN.size()*nndof)
{
std::cerr <<" *** ASMbase::injectNodeVec:: Invalid patch vector, size = "
<< nodeVec.size() <<" < "<< MLGN.size()*nndof << std::endl;
<< nodeVec.size() <<" != "<< MLGN.size()*nndof << std::endl;
return false;
}

53
src/SIM/SIMbase.C
View File

@ -1036,17 +1036,19 @@ bool SIMbase::writeGlvV (const Vector& vec, const char* fieldName,
bool SIMbase::writeGlvS (const Vector& psol, const int* nViz,
int iStep, int& nBlock, double time,
bool psolOnly, const char* vecName)
char psolOnly, const char* pvecName,
int idBlock, int psolComps)
{
if (psol.empty())
return true;
else if (!myVtf)
return false;
if (!psolOnly)
if (!psolOnly && myProblem)
myProblem->initResultPoints(time);
Matrix field;
Vector lovec;
size_t i, j;
int geomID = 0;
const size_t pMAX = 6;
@ -1070,11 +1072,14 @@ bool SIMbase::writeGlvS (const Vector& psol, const int* nViz,
// 1. Evaluate primary solution variables
myModel[i]->extractNodeVec(psol,myProblem->getSolution());
if (!myModel[i]->evalSolution(field,myProblem->getSolution(),nViz))
Vector& locvec = myProblem ? myProblem->getSolution() : lovec;
myModel[i]->extractNodeVec(psol,locvec,psolComps);
if (!myModel[i]->evalSolution(field,locvec,nViz))
return false;
if (!myVtf->writeVres(field,++nBlock,++geomID,nVcomp))
if (psolOnly > 1)
geomID++;
else if (!myVtf->writeVres(field,++nBlock,++geomID,nVcomp))
return false;
else
vID[0].push_back(nBlock);
@ -1085,7 +1090,7 @@ bool SIMbase::writeGlvS (const Vector& psol, const int* nViz,
else
dID[j].push_back(nBlock);
if (psolOnly) continue; // skip secondary solution
if (psolOnly || !myProblem) continue; // skip secondary solution
// 2. Direct evaluation of secondary solution variables
@ -1154,22 +1159,32 @@ bool SIMbase::writeGlvS (const Vector& psol, const int* nViz,
// Write result block identifications
bool ok = true;
int idBlock = 10;
if (vecName)
ok = myVtf->writeVblk(vID[0],vecName,idBlock,iStep);
else
ok = myVtf->writeDblk(vID[0],"Solution",idBlock,iStep);
if (!ok) return false;
if (!vID[0].empty())
if (pvecName)
ok = myVtf->writeVblk(vID[0],pvecName,idBlock,iStep);
else
ok = myVtf->writeDblk(vID[0],"Solution",idBlock,iStep);
if (scalarEq && !psolOnly)
if (!myVtf->writeVblk(vID[1],"Flux",idBlock+1,iStep))
return false;
if (ok && !vID[1].empty())
ok = myVtf->writeVblk(vID[1],"Flux",idBlock+1,iStep);
for (j = 0; j < pMAX && !dID[j].empty(); j++)
if (!myVtf->writeSblk(dID[j],myProblem->getField1Name(j),++idBlock,iStep))
return false;
std::string pname;
if (!myProblem)
{
pname = pvecName ? pvecName : "Solution";
if (psolOnly < 2) pname += "_w";
}
if (psolOnly) return true;
for (j = 0; j < pMAX && !dID[j].empty() && ok; j++)
{
if (myProblem)
pname = myProblem->getField1Name(j);
else if (psolOnly < 2)
(*pname.rbegin()) ++;
ok = myVtf->writeSblk(dID[j],pname.c_str(),++idBlock,iStep);
}
if (psolOnly || !myProblem || !ok) return ok;
size_t nf = myProblem->getNoFields(2);
for (i = j = 0; i < nf && !sID[j].empty(); i++, j++)

11
src/SIM/SIMbase.h
View File

@ -296,15 +296,18 @@ public:
//! \param nBlock Running result block counter
//! \param[in] time Load/time step parameter
//! \param[in] psolOnly If \e true, skip secondary solution field evaluation
//! \param[in] vecName Optional name of the primary vector field solution
//! \param[in] pvecName Optional name of the primary vector field solution
//! \param[in] idBlock Starting value of result block numbering
//! \param[in] psolComps Optional number of primary solution components
//!
//! \details The primary vector field solution is written as a deformation
//! plot if \a vecName is NULL. Otherwise, it is written as a named vector
//! plot if \a pvecName is NULL. Otherwise, it is written as a named vector
//! field. If an analytical solution is provided, the exact secondary solution
//! fields are written to the VTF-file as well.
virtual bool writeGlvS(const Vector& psol, const int* nViz, int iStep,
int& nBlock, double time = 0.0, bool psolOnly = false,
const char* vecName = 0);
int& nBlock, double time = 0.0, char psolOnly = 0,
const char* pvecName = 0, int idBlock = 10,
int psolComps = 0);
//! \brief Writes a mode shape and associated eigenvalue to the VTF-file.
//! \details The eigenvalue is used as a label on the step state info

76
src/Utility/DataExporter.C
View File

@ -3,15 +3,26 @@
#include "DataExporter.h"
#include <iostream>
#include <algorithm>
#ifdef PARALLEL_PETSC
#include <mpi.h>
#endif
DataExporter::DataExporter(bool dynWrts) : deleteWriters(dynWrts), m_level(-1)
DataWriter::DataWriter (const std::string& name) : m_name(name)
{
#ifdef PARALLEL_PETSC
MPI_Comm_size(MPI_COMM_WORLD,&m_size);
MPI_Comm_rank(MPI_COMM_WORLD,&m_rank);
#else
m_size = 1;
m_rank = 0;
#endif
}
DataExporter::~DataExporter()
DataExporter::~DataExporter ()
{
if (deleteWriters)
if (m_delete)
for (size_t i = 0; i < m_writers.size(); i++)
delete m_writers[i];
}
@ -40,6 +51,7 @@ bool DataExporter::registerWriter(DataWriter* writer)
return true;
}
bool DataExporter::setFieldValue(const std::string& name,
void* data, void* data2)
{
@ -52,10 +64,11 @@ bool DataExporter::setFieldValue(const std::string& name,
return true;
}
bool DataExporter::dumpTimeLevel()
{
if (m_level == -1)
m_level = getWritersTimeLevel()+1;
m_level = this->getWritersTimeLevel()+1;
std::map<std::string,FileEntry>::iterator it;
std::vector<DataWriter*>::iterator it2;
@ -72,7 +85,7 @@ bool DataExporter::dumpTimeLevel()
(*it2)->writeSIM(m_level,*it);
break;
default:
std::cout <<"DataExporter: Invalid field type registered, skipping"
std::cerr <<"DataExporter: Invalid field type registered, skipping"
<< std::endl;
break;
}
@ -84,52 +97,59 @@ bool DataExporter::dumpTimeLevel()
return true;
}
bool DataExporter::loadTimeLevel(int level, DataWriter* input)
{
if (!input && m_writers.empty())
return false;
if (!input)
input = m_writers.front();
if (level == -1)
level = m_level = input->getLastTimeLevel();
if (level == -1)
return false;
bool DataExporter::loadTimeLevel (int level, DataWriter* input)
{
if (!input)
if (m_writers.empty())
return false;
else
input = m_writers.front();
if (level == -1)
if ((m_level = input->getLastTimeLevel()) < 0)
return false;
else
level = m_level;
bool ok = true;
input->openFile(level);
std::map<std::string,FileEntry>::iterator it;
for (it = m_entry.begin(); it != m_entry.end(); ++it)
{
for (it = m_entry.begin(); it != m_entry.end() && ok; ++it)
if (!it->second.data)
return false;
switch (it->second.field)
ok = false;
else switch (it->second.field)
{
case VECTOR:
input->readVector(level,*it);
ok = input->readVector(level,*it);
break;
case SIM:
input->readSIM(level,*it);
ok = input->readSIM(level,*it);
break;
default:
std::cout << "DataExporter::loadTimeLevel: Invalid field type "
<< "registered, skipping" << std::endl;
ok = false;
std::cerr <<" *** DataExporter: Invalid field type registered "
<< it->second.field << std::endl;
break;
}
}
input->closeFile(level);
m_level++;
return true;
return ok;
}
int DataExporter::getTimeLevel()
int DataExporter::getTimeLevel ()
{
if (m_level == -1)
m_level = getWritersTimeLevel();
m_level = this->getWritersTimeLevel();
return m_level;
}
int DataExporter::getWritersTimeLevel()
int DataExporter::getWritersTimeLevel () const
{
std::vector<int> levels;
std::vector<DataWriter*>::const_iterator it2;

88
src/Utility/DataExporter.h
View File

@ -9,51 +9,55 @@
class DataWriter;
class DataExporter {
public:
enum FieldType {
VECTOR,
SIM
};
class DataExporter
{
public:
enum FieldType {
VECTOR,
SIM
};
struct FileEntry {
std::string description;
FieldType field;
int size;
void* data;
void* data2;
};
struct FileEntry {
std::string description;
FieldType field;
int size;
void* data;
void* data2;
};
DataExporter(bool dynamicWriters = false);
virtual ~DataExporter();
//! \brief Default constructor.
//! \param[in] dynWriters If \e true, delete the writers on destruction.
DataExporter(bool dynWriters = false) : m_delete(dynWriters), m_level(-1) {}
//! \brief The desctructor deletes the writers if \a dynWriters was \e true.
~DataExporter();
//! \brief Registers an entry for storage.
//! param[in] name Name of entry
//! param[in] description Description of entry
//! param[in] Type of entry
//! param[in] size set to number of entries in an array,
// the time level to use for SIM
bool registerField(const std::string& name,
const std::string& description,
FieldType field, size_t size=0);
//! \brief Registers an entry for storage.
//! param[in] name Name of entry
//! param[in] description Description of entry
//! param[in] field Type of entry
//! param[in] size Number of entries in an array,
//! the time level to use for SIM
bool registerField(const std::string& name,
const std::string& description,
FieldType field, size_t size=0);
bool registerWriter(DataWriter* writer);
bool registerWriter(DataWriter* writer);
bool setFieldValue(const std::string& name, void* data, void* data2=NULL);
bool setFieldValue(const std::string& name, void* data, void* data2=NULL);
bool dumpTimeLevel();
bool dumpTimeLevel();
//! \brief Loads last time level with first registered writer by default.
bool loadTimeLevel(int level=-1, DataWriter* input=NULL);
int getTimeLevel();
//! \brief Loads last time level with first registered writer by default.
bool loadTimeLevel(int level=-1, DataWriter* input=NULL);
int getTimeLevel();
protected:
int getWritersTimeLevel();
protected:
int getWritersTimeLevel() const;
std::map<std::string,FileEntry> m_entry;
std::vector<DataWriter*> m_writers;
bool deleteWriters;
int m_level;
std::map<std::string,FileEntry> m_entry;
std::vector<DataWriter*> m_writers;
bool m_delete;
int m_level;
};
@ -62,7 +66,7 @@ typedef std::pair<std::string,DataExporter::FileEntry> DataEntry;
class DataWriter
{
protected:
DataWriter() {}
DataWriter(const std::string& name);
public:
virtual ~DataWriter() {}
@ -73,7 +77,13 @@ public:
virtual void closeFile(int level) = 0;
virtual void writeVector(int level, const DataEntry& entry) = 0;
virtual void readVector(int level, const DataEntry& entry) = 0;
virtual bool readVector(int level, const DataEntry& entry) = 0;
virtual void writeSIM(int level, const DataEntry& entry) = 0;
virtual void readSIM(int level, const DataEntry& entry) = 0;
virtual bool readSIM(int level, const DataEntry& entry) = 0;
protected:
std::string m_name; //!< File name
int m_size; //!< Number of MPI nodes (processors)
int m_rank; //!< MPI rank (processor ID)
};

108
src/Utility/HDF5Writer.C
View File

@ -4,7 +4,6 @@
#include "SIMbase.h"
#include "IntegrandBase.h"
#include <sstream>
#include <numeric>
#include <sys/stat.h>
#ifdef HAS_HDF5
@ -16,30 +15,23 @@
#endif
HDF5Writer::HDF5Writer (const std::string& name) : m_hdf5(name+".hdf5")
HDF5Writer::HDF5Writer (const std::string& name, bool append)
: DataWriter(name+".hdf5")
{
#ifdef HAS_HDF5
struct stat temp;
// file already exists - open and find the next group
if (stat(m_hdf5.c_str(),&temp) == 0)
if (append && stat(m_name.c_str(),&temp) == 0)
m_flag = H5F_ACC_RDWR;
else
m_flag = H5F_ACC_TRUNC;
#endif
#ifdef PARALLEL_PETSC
MPI_Comm_rank(MPI_COMM_WORLD,&m_rank);
MPI_Comm_size(MPI_COMM_WORLD,&m_size);
#else
m_rank = 0;
#endif
}
HDF5Writer::~HDF5Writer()
{
}
int HDF5Writer::getLastTimeLevel()
{
int result = 0;
#ifdef HAS_HDF5
if (m_flag == H5F_ACC_TRUNC)
return -1;
@ -47,12 +39,11 @@ int HDF5Writer::getLastTimeLevel()
hid_t acc_tpl = H5P_DEFAULT;
#ifdef PARALLEL_PETSC
MPI_Info info = MPI_INFO_NULL;
acc_tpl = H5Pcreate (H5P_FILE_ACCESS);
acc_tpl = H5Pcreate(H5P_FILE_ACCESS);
H5Pset_fapl_mpio(acc_tpl, MPI_COMM_WORLD, info);
#endif
m_file = H5Fopen(m_hdf5.c_str(),m_flag,acc_tpl);
int result=0;
m_file = H5Fopen(m_name.c_str(),m_flag,acc_tpl);
while (1) {
std::stringstream str;
str << '/' << result;
@ -61,11 +52,8 @@ int HDF5Writer::getLastTimeLevel()
result++;
}
H5Fclose(m_file);
return result-1;
#else
return -1;
#endif
return result-1;
}
void HDF5Writer::openFile(int level)
@ -78,9 +66,9 @@ void HDF5Writer::openFile(int level)
H5Pset_fapl_mpio(acc_tpl, MPI_COMM_WORLD, info);
#endif
if (m_flag == H5F_ACC_TRUNC)
m_file = H5Fcreate(m_hdf5.c_str(),m_flag,H5P_DEFAULT,acc_tpl);
m_file = H5Fcreate(m_name.c_str(),m_flag,H5P_DEFAULT,acc_tpl);
else
m_file = H5Fopen(m_hdf5.c_str(),m_flag,acc_tpl);
m_file = H5Fopen(m_name.c_str(),m_flag,acc_tpl);
std::stringstream str;
str << '/' << level;
@ -153,9 +141,10 @@ void HDF5Writer::writeArray(int group, const std::string& name,
#endif
}
void HDF5Writer::readVector(int level, const DataEntry& entry)
bool HDF5Writer::readVector(int level, const DataEntry& entry)
{
// readArray(level,entry.first,entry.second.size,entry.second.data);
return true;
}
void HDF5Writer::writeVector(int level, const DataEntry& entry)
@ -163,14 +152,15 @@ void HDF5Writer::writeVector(int level, const DataEntry& entry)
writeArray(level,entry.first,entry.second.size,entry.second.data);
}
void HDF5Writer::readSIM(int level, const DataEntry& entry)
bool HDF5Writer::readSIM(int level, const DataEntry& entry)
{
bool ok = true;
#ifdef HAS_HDF5
SIMbase* sim = static_cast<SIMbase*>(entry.second.data);
Vector* sol = static_cast<Vector*>(entry.second.data2);
if (!sol) return;
if (!sol) return false;
for (int i = 0; i < sim->getNoPatches(); ++i) {
for (int i = 0; i < sim->getNoPatches() && ok; ++i) {
std::stringstream str;
str << level;
str << '/';
@ -180,12 +170,42 @@ void HDF5Writer::readSIM(int level, const DataEntry& entry)
if (loc > 0) {
double* tmp = NULL; int siz = 0;
readArray(group2,entry.first,siz,tmp);
sim->injectPatchSolution(*sol,Vector(tmp,siz),loc-1,entry.second.size);
ok = sim->injectPatchSolution(*sol,Vector(tmp,siz),
loc-1,entry.second.size);
delete[] tmp;
}
H5Gclose(group2);
}
#endif
return ok;
}
bool HDF5Writer::readField(int level, const std::string& name,
Vector& vec, SIMbase* sim, int components)
{
bool ok = true;
openFile(level);
#ifdef HAS_HDF5
vec.resize(sim->getNoNodes()*components);
for (int i = 0; i < sim->getNoPatches() && ok; ++i) {
std::stringstream str;
str << level;
str << '/';
str << i+1;
hid_t group2 = H5Gopen2(m_file,str.str().c_str(),H5P_DEFAULT);
int loc = sim->getLocalPatchIndex(i+1);
if (loc > 0) {
double* tmp = NULL; int siz = 0;
readArray(group2,name,siz,tmp);
ok = sim->injectPatchSolution(vec,Vector(tmp,siz),
loc-1,components);
delete[] tmp;
}
H5Gclose(group2);
}
#endif
closeFile(level);
return ok;
}
void HDF5Writer::writeSIM(int level, const DataEntry& entry)
@ -208,17 +228,29 @@ void HDF5Writer::writeSIM(int level, const DataEntry& entry)
group2 = H5Gcreate2(m_file,str.str().c_str(),0,H5P_DEFAULT,H5P_DEFAULT);
int loc = sim->getLocalPatchIndex(i+1);
if (loc > 0) // we own the patch
{
sim->extractPatchSolution(*sol,loc-1);
if (entry.second.size == -1) {
Matrix field;
if (loc > 0)
sim->evalSecondarySolution(field,loc-1);
for (size_t j = 0; j < prob->getNoFields(2); ++j)
writeArray(group2,prob->getField2Name(j),field.cols(),
field.getRow(j+1).ptr());
Vector& psol = const_cast<Integrand*>(prob)->getSolution();
writeArray(group2,entry.first,psol.size(),psol.ptr());
// TODO: For mixed methods we need to output each primary field
// TODO: The above is sufficient for simulation restarts.
// TODO: For mixed methods we need to output each primary field
// TODO: separately in addition, with reference to correct spline basis.
if (entry.second.size == -1) {
Matrix field;
sim->evalSecondarySolution(field,loc-1);
for (size_t j = 0; j < prob->getNoFields(2); ++j)
writeArray(group2,prob->getField2Name(j),field.cols(),
field.getRow(j+1).ptr());
}
}
else // must write empty dummy records for the other patches
{
double dummy;
writeArray(group2,entry.first,0,&dummy);
for (size_t j = 0; j < prob->getNoFields(2); ++j)
writeArray(group2,prob->getField2Name(j),0,&dummy);
}
Vector& psol = const_cast<Integrand*>(prob)->getSolution();
writeArray(group2, entry.first, loc > 0 ? psol.size() : 0, psol.ptr());
H5Gclose(group2);
}
#endif
@ -226,15 +258,15 @@ void HDF5Writer::writeSIM(int level, const DataEntry& entry)
bool HDF5Writer::checkGroupExistence(int parent, const char* path)
{
bool result = false;
#ifdef HAS_HDF5
// turn off errors to avoid cout spew
herr_t (*old_func)(void*);
void* old_client_data;
H5Eget_auto(&old_func,&old_client_data);
H5Eset_auto(NULL,NULL);
bool result =H5Gget_objinfo((hid_t)parent,path,0,NULL) == 0;
result = H5Gget_objinfo((hid_t)parent,path,0,NULL) == 0;
H5Eset_auto(old_func,old_client_data);
return result;
#endif
return false;
return result;
}

48
src/Utility/HDF5Writer.h
View File

@ -3,35 +3,37 @@
#pragma once
#include "DataExporter.h"
#include "MatVec.h"
class SIMbase;
class HDF5Writer : public DataWriter {
public:
class HDF5Writer : public DataWriter
{
public:
HDF5Writer(const std::string& name, bool append = false);
virtual ~HDF5Writer() {}
HDF5Writer(const std::string& name);
virtual ~HDF5Writer();
virtual int getLastTimeLevel();
int getLastTimeLevel();
virtual void openFile(int level);
virtual void closeFile(int level);
void openFile(int level);
void closeFile(int level);
virtual void writeVector(int level, const DataEntry& entry);
virtual bool readVector(int level, const DataEntry& entry);
virtual void writeSIM(int level, const DataEntry& entry);
virtual bool readSIM(int level, const DataEntry& entry);
void writeVector(int level, const DataEntry& entry);
void readVector(int level, const DataEntry& entry);
void writeSIM(int level, const DataEntry& entry);
void readSIM(int level, const DataEntry& entry);
bool readField(int level, const std::string& name,
Vector& vec, SIMbase* sim, int components);
protected:
void writeArray(int group, const std::string& name,
int len, void* data);
void readArray(int group, const std::string& name,
int& len, double*& data);
bool checkGroupExistence(int parent, const char* group);
protected:
void writeArray(int group, const std::string& name,
int len, void* data);
void readArray(int group, const std::string& name,
int& len, double*& data);
bool checkGroupExistence(int parent, const char* group);
std::string m_hdf5;
int m_file;
unsigned int m_flag;
int m_rank; // MPI rank
int m_size; // number of MPI nodes
int m_file;
unsigned int m_flag;
};

6
src/Utility/Utilities.h
View File

@ -23,10 +23,10 @@
namespace utl
{
//! \brief Parses a character string into an integer or an integer range.
//! \param values The integer value(s) is(are) appended to this vector
//! \param argv Character string with integer data
//! \param values The integer value(s) is/are appended to this vector
//! \param[in] argv Character string with integer data
//!
//! \details An integer range is recognised through the syntax <i>:<j>.
//! \details An integer range is recognised through the syntax \a i:j.
void parseIntegers(std::vector<int>& values, const char* argv);
//! \brief Reads one line, ignoring comment lines and leading blanks.

77
src/Utility/XMLWriter.C
View File

@ -3,31 +3,23 @@
#include "XMLWriter.h"
#include "SIMbase.h"
#include "IntegrandBase.h"
#ifdef PARALLEL_PETSC
#include <mpi.h>
#endif
#include "StringUtils.h"
#include "tinyxml.h"
#include <fstream>
#include <cstdio>
XMLWriter::XMLWriter(const std::string& name) : m_xml(name+".xml")
XMLWriter::XMLWriter(const std::string& name) : DataWriter(name+".xml")
{
m_doc = NULL;
#ifdef PARALLEL_PETSC
MPI_Comm_rank(MPI_COMM_WORLD,&m_rank);
MPI_Comm_size(MPI_COMM_WORLD,&m_size);
#else
m_rank = 0;
#endif
m_node = NULL;
}
XMLWriter::~XMLWriter()
{
}
int XMLWriter::getLastTimeLevel()
{
int result=-1;
TiXmlDocument doc(m_xml.c_str());
TiXmlDocument doc(m_name.c_str());
doc.LoadFile();
TiXmlHandle handle(&doc);
TiXmlElement* levels = handle.FirstChild("info").
@ -38,10 +30,12 @@ int XMLWriter::getLastTimeLevel()
return result;
}
void XMLWriter::openFile(int level)
{
if (m_rank != 0)
return;
m_doc = new TiXmlDocument;
TiXmlElement element("info");
m_node = m_doc->InsertEndChild(element);
@ -59,13 +53,37 @@ void XMLWriter::closeFile(int level)
TiXmlText value(temp);
pNewNode->InsertEndChild(value);
m_doc->SaveFile(m_xml);
m_doc->SaveFile(m_name);
delete m_doc;
m_doc = NULL;
}
void XMLWriter::readVector(int level, const DataEntry& entry)
void XMLWriter::readInfo()
{
TiXmlDocument doc(m_name);
doc.LoadFile();
TiXmlHandle handle(&doc);
TiXmlElement* elem = handle.FirstChild("info").
FirstChild("entry").ToElement();
while (elem) {
if (strcasecmp(elem->Attribute("type"),"field") == 0) {
Entry entry;
entry.name = elem->Attribute("name");
entry.description = elem->Attribute("description");
entry.patches = atoi(elem->Attribute("patches"));
entry.components = atoi(elem->Attribute("components"));
entry.patchfile = elem->Attribute("patchfile");
m_entry.push_back(entry);
}
elem = elem->NextSiblingElement("entry");
}
}
bool XMLWriter::readVector(int level, const DataEntry& entry)
{
return true;
}
void XMLWriter::writeVector(int level, const DataEntry& entry)
@ -81,8 +99,10 @@ void XMLWriter::writeVector(int level, const DataEntry& entry)
m_node->InsertEndChild(element);
}
void XMLWriter::readSIM(int level, const DataEntry& entry)
bool XMLWriter::readSIM(int level, const DataEntry& entry)
{
return true;
}
void XMLWriter::writeSIM(int level, const DataEntry& entry)
@ -90,25 +110,38 @@ void XMLWriter::writeSIM(int level, const DataEntry& entry)
if (m_rank != 0)
return;
// Assume that all fields use the same basis as the geometry.
// TODO: Not true for mixed methods, fix later...
std::string g2file(m_name);
std::ofstream os(replaceAll(g2file,".xml",".g2").c_str());
static_cast<SIMbase*>(entry.second.data)->dumpGeometry(os);
int nPatch = static_cast<SIMbase*>(entry.second.data)->getNoPatches();
const Integrand* prb = static_cast<SIMbase*>(entry.second.data)->getProblem();
// primary solution
addField(entry.first,entry.second.description,"field",nPatch);
addField(entry.first,entry.second.description,"field",g2file,
prb->getNoFields(1),nPatch);
// secondary solutions
if (entry.second.size == -1)
for (size_t j = 0; j < prb->getNoFields(2); j++)
addField(prb->getField2Name(j),"secondary","field",nPatch);
addField(prb->getField2Name(j),"secondary","field",g2file,1,nPatch);
}
void XMLWriter::addField(const std::string& name, const std::string& description,
const std::string& type, int patches)
void XMLWriter::addField (const std::string& name,
const std::string& description,
const std::string& type,
const std::string& patchfile,
int components, int patches)
{
TiXmlElement element("entry");
element.SetAttribute("name",name.c_str());
element.SetAttribute("description",description.c_str());
element.SetAttribute("type",type.c_str());
element.SetAttribute("patches",patches);
element.SetAttribute("components",components);
element.SetAttribute("patchfile",patchfile.c_str());
m_node->InsertEndChild(element);
}

53
src/Utility/XMLWriter.h
View File

@ -3,32 +3,45 @@
#pragma once
#include "DataExporter.h"
#include "tinyxml.h"
class TiXmlDocument;
class TiXmlNode;
class XMLWriter : public DataWriter {
public:
XMLWriter(const std::string& name);
virtual ~XMLWriter();
class XMLWriter : public DataWriter
{
public:
struct Entry {
std::string name;
std::string description;
std::string patchfile;
int patches;
int components;
};
int getLastTimeLevel();
XMLWriter(const std::string& name);
virtual ~XMLWriter() {}
void openFile(int level);
void closeFile(int level);
virtual int getLastTimeLevel();
void writeVector(int level, const DataEntry& entry);
void readVector(int level, const DataEntry& entry);
void writeSIM(int level, const DataEntry& entry);
void readSIM(int level, const DataEntry& entry);
void readInfo();
const std::vector<Entry>& getEntries() const { return m_entry; }
protected:
void addField(const std::string& name, const std::string& description,
const std::string& type, int patches);
std::string m_xml;
virtual void openFile(int level);
virtual void closeFile(int level);
int m_rank; // MPI rank
int m_size; // number of MPI nodes
virtual void writeVector(int level, const DataEntry& entry);
virtual bool readVector(int level, const DataEntry& entry);
virtual void writeSIM(int level, const DataEntry& entry);
virtual bool readSIM(int level, const DataEntry& entry);
TiXmlDocument* m_doc;
TiXmlNode* m_node;
protected:
void addField(const std::string& name, const std::string& description,
const std::string& type, const std::string& patchfile,
int components, int patches);
std::vector<Entry> m_entry;
TiXmlDocument* m_doc;
TiXmlNode* m_node;
};