Files
LBPM/visit/avtLBMFileFormat.C

485 lines
18 KiB
C

/*****************************************************************************
*
* Copyright (c) 2000 - 2013, Lawrence Livermore National Security, LLC
* Produced at the Lawrence Livermore National Laboratory
* LLNL-CODE-442911
* All rights reserved.
*
* This file is part of VisIt. For details, see https://visit.llnl.gov/. The
* full copyright notice is contained in the file COPYRIGHT located at the root
* of the VisIt distribution or at http://www.llnl.gov/visit/copyright.html.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* - Redistributions of source code must retain the above copyright notice,
* this list of conditions and the disclaimer below.
* - Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the disclaimer (as noted below) in the
* documentation and/or other materials provided with the distribution.
* - Neither the name of the LLNS/LLNL nor the names of its contributors may
* be used to endorse or promote products derived from this software without
* specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL LAWRENCE LIVERMORE NATIONAL SECURITY,
* LLC, THE U.S. DEPARTMENT OF ENERGY OR CONTRIBUTORS BE LIABLE FOR ANY
* DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
* SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
* CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
* DAMAGE.
*
*****************************************************************************/
// ************************************************************************* //
// avtLBMFileFormat.C //
// ************************************************************************* //
#include <avtLBMFileFormat.h>
#include <DebugStream.h>
#include <string>
// LBPM headers
#include "IO/Reader.h"
#include "IO/IOHelpers.h"
#include "IO/PIO.h"
#include "common/Utilities.h"
// vtk headers
#include <vtkFloatArray.h>
#include <vtkRectilinearGrid.h>
#include <vtkStructuredGrid.h>
#include <vtkUnstructuredGrid.h>
#include <vtkPointSet.h>
#include <vtkCellType.h>
#include <vtkPolyData.h>
#include <vtkCellArray.h>
#include <avtDatabaseMetaData.h>
#include <vtkHelpers.h>
#include <DBOptionsAttributes.h>
#include <Expression.h>
#include <InvalidVariableException.h>
#if defined(WIN32) || defined(_WIN32) || defined(WIN64) || defined(_WIN64)
#define USE_WINDOWS
#elif defined(__APPLE__)
#define USE_MAC
#else
#define USE_LINUX
#endif
#ifdef USE_WINDOWS
#include <windows.h>
#include <stdio.h>
#define TIME_TYPE LARGE_INTEGER
#define get_time(x) QueryPerformanceCounter(x)
#define get_frequency(f) QueryPerformanceFrequency(f)
#define get_diff(start,end,f) \
static_cast<double>(end.QuadPart-start.QuadPart)/static_cast<double>(f.QuadPart)
#elif defined(USE_LINUX) || defined(USE_MAX)
#include <sys/time.h>
#define TIME_TYPE timeval
#define get_time(x) gettimeofday(x,NULL);
#define get_frequency(f) (*f=timeval())
#define get_diff(start,end,f) 1e-6*static_cast<double>( \
0xF4240*(static_cast<int64_t>(end.tv_sec)-static_cast<int64_t>(start.tv_sec)) + \
(static_cast<int64_t>(end.tv_usec)-static_cast<int64_t>(start.tv_usec)) )
#else
#error Unknown OS
#endif
// Output streams
static IO::ParallelStreamBuffer DebugStreamBuffer1;
static IO::ParallelStreamBuffer DebugStreamBuffer2;
std::ostream DebugStream1(&DebugStreamBuffer1);
std::ostream DebugStream2(&DebugStreamBuffer2);
// ****************************************************************************
// Method: avtLBMFileFormat constructor
//
// Programmer: mbt -- generated by xml2avt
// Creation: Wed Jul 9 10:59:08 PDT 2014
//
// ****************************************************************************
avtLBMFileFormat::avtLBMFileFormat(const char *filename)
: avtMTMDFileFormat(filename)
{
// Set abort behavior
Utilities::setAbortBehavior(true,true,true);
Utilities::setErrorHandlers();
// Set debug streams
DebugStreamBuffer1.setOutputStream( &DebugStream::Stream1() );
DebugStreamBuffer2.setOutputStream( &DebugStream::Stream2() );
DebugStreamBuffer1.setOutputStream( &std::cout );
// Get the path to the input file
std::string file(filename);
size_t k1 = file.rfind(47);
size_t k2 = file.rfind(92);
if ( k1==std::string::npos ) { k1=0; }
if ( k2==std::string::npos ) { k2=0; }
d_path = file.substr(0,std::max(k1,k2));
// Load the summary file
DebugStream1 << "Loading " << filename << std::endl;
d_timesteps = IO::readTimesteps(filename);
// Read the mesh dabases
d_database.clear();
d_database.resize(d_timesteps.size());
for (size_t i=0; i<d_timesteps.size(); i++)
d_database[i] = IO::getMeshList(d_path,d_timesteps[i]);
}
// ****************************************************************************
// Method: avtEMSTDFileFormat::GetNTimesteps
//
// Purpose:
// Tells the rest of the code how many timesteps there are in this file.
//
// Programmer: mbt -- generated by xml2avt
// Creation: Wed Jul 9 10:59:08 PDT 2014
//
// ****************************************************************************
int
avtLBMFileFormat::GetNTimesteps(void)
{
DebugStream2 << "avtLBMFileFormat::GetNTimesteps" << std::endl;
return static_cast<int>(d_timesteps.size());
}
// ****************************************************************************
// Method: avtLBMFileFormat::FreeUpResources
//
// Purpose:
// When VisIt is done focusing on a particular timestep, it asks that
// timestep to free up any resources (memory, file descriptors) that
// it has associated with it. This method is the mechanism for doing
// that.
//
// Programmer: mbt -- generated by xml2avt
// Creation: Wed Jul 9 10:59:08 PDT 2014
//
// ****************************************************************************
void
avtLBMFileFormat::FreeUpResources(void)
{
DebugStream1 << "avtLBMFileFormat::FreeUpResources" << std::endl;
std::map<std::string,vtkObjectBase*>::iterator it;
for ( it=d_meshcache.begin(); it!=d_meshcache.end(); ++it ) {
DebugStream2 << " deleting: " << it->first << std::endl;
vtkObjectBase* obj = it->second;
it->second = NULL;
if ( obj!=NULL )
obj->Delete();
}
DebugStream2 << " finished" << std::endl;
}
// ****************************************************************************
// Method: avtLBMFileFormat::PopulateDatabaseMetaData
//
// Purpose:
// This database meta-data object is like a table of contents for the
// file. By populating it, you are telling the rest of VisIt what
// information it can request from you.
//
// Programmer: mbt -- generated by xml2avt
// Creation: Wed Jul 9 10:59:08 PDT 2014
//
// ****************************************************************************
void
avtLBMFileFormat::PopulateDatabaseMetaData(avtDatabaseMetaData *md, int timestate)
{
DebugStream1 << "avtLBMFileFormat::PopulateDatabaseMetaData: " << timestate << std::endl;
// Add the mesh domains to the meta data
const std::vector<IO::MeshDatabase> database = d_database[timestate];
for (size_t i=0; i<database.size(); i++) {
DebugStream2 << " Adding " << database[i].name << std::endl;
avtMeshMetaData *mmd = new avtMeshMetaData;
mmd->name = database[i].name;
mmd->meshType = vtkMeshType(database[i].type);
mmd->spatialDimension = 3;
mmd->topologicalDimension = vtkTopDim(database[i].type);
if ( mmd->meshType==AVT_SURFACE_MESH )
mmd->topologicalDimension = 2;
mmd->numBlocks = database[i].domains.size();
mmd->blockNames.resize( mmd->numBlocks );
for (int j=0; j<mmd->numBlocks; j++)
mmd->blockNames[j] = database[i].domains[j].name;
md->Add(mmd);
// Add expressions for the coordinates
const char *xyz[3] = {"x","y","z"};
for(int j=0; j<3; j++) {
Expression expr;
char expdef[100], expname[100];
sprintf(expdef,"coord(<%s>)[%i]",mmd->name.c_str(),j);
sprintf(expname,"%s/%s",xyz[j],mmd->name.c_str());
expr.SetName(expname);
expr.SetDefinition(expdef);
md->AddExpression(&expr);
}
// Add the variables
for (size_t j=0; j<database[i].variables.size(); j++) {
IO::VariableDatabase variable = database[i].variables[j];
std::string varname = variable.name + "/" + mmd->name;
avtCentering center = AVT_UNKNOWN_CENT;
if ( variable.type==IO::NodeVariable ) {
center = AVT_NODECENT;
} else if ( variable.type==IO::SurfaceVariable ) {
center = AVT_ZONECENT;
} else if ( variable.type==IO::VolumeVariable ) {
center = AVT_ZONECENT;
}
if ( variable.dim==1 ) {
AddScalarVarToMetaData( md, varname, mmd->name, center );
} else if ( variable.dim==3 ) {
AddVectorVarToMetaData( md, varname, mmd->name, center, variable.dim );
} else if ( variable.dim==9 ) {
AddTensorVarToMetaData( md, varname, mmd->name, center, variable.dim );
}
}
}
DebugStream2 << " Finished" << std::endl;
//
// CODE TO ADD A MATERIAL
//
// string mesh_for_mat = meshname; // ??? -- could be multiple meshes
// string matname = ...
// int nmats = ...;
// vector<string> mnames;
// for (int i = 0 ; i < nmats ; i++)
// {
// char str[32];
// sprintf(str, "mat%d", i);
// -- or --
// strcpy(str, "Aluminum");
// mnames.push_back(str);
// }
//
// Here's the call that tells the meta-data object that we have a mat:
//
// AddMaterialToMetaData(md, matname, mesh_for_mat, nmats, mnames);
//
//
// Here's the way to add expressions:
//Expression momentum_expr;
//momentum_expr.SetName("momentum");
//momentum_expr.SetDefinition("{u, v}");
//momentum_expr.SetType(Expression::VectorMeshVar);
//md->AddExpression(&momentum_expr);
//Expression KineticEnergy_expr;
//KineticEnergy_expr.SetName("KineticEnergy");
//KineticEnergy_expr.SetDefinition("0.5*(momentum*momentum)/(rho*rho)");
//KineticEnergy_expr.SetType(Expression::ScalarMeshVar);
//md->AddExpression(&KineticEnergy_expr);
//
}
// ****************************************************************************
// Method: avtLBMFileFormat::GetMesh
//
// Purpose:
// Gets the mesh associated with this file. The mesh is returned as a
// derived type of vtkDataSet (ie vtkRectilinearGrid, vtkStructuredGrid,
// vtkUnstructuredGrid, etc).
//
// Arguments:
// timestate The index of the timestate. If GetNTimesteps returned
// 'N' time steps, this is guaranteed to be between 0 and N-1.
// domain The index of the domain. If there are NDomains, this
// value is guaranteed to be between 0 and NDomains-1,
// regardless of block origin.
// meshname The name of the mesh of interest. This can be ignored if
// there is only one mesh.
//
// Programmer: mbt -- generated by xml2avt
// Creation: Wed Jul 9 10:59:08 PDT 2014
//
// ****************************************************************************
vtkDataSet *
avtLBMFileFormat::GetMesh(int timestate, int domain, const char *meshname)
{
DebugStream1 << "avtLBMFileFormat::GetMesh - " << meshname
<< "," << timestate << "," << domain << std::endl;
TIME_TYPE start, stop, freq;
get_frequency(&freq);
get_time(&start);
// Check if we have a cached copy of the mesh
char cache_name[1000];
sprintf(cache_name,"%i-%i-%s",timestate,domain,meshname);
if ( d_meshcache.find(cache_name)!=d_meshcache.end() )
return dynamic_cast<vtkDataSet*>(d_meshcache.find(cache_name)->second);
// Read the mesh
std::shared_ptr<IO::Mesh> mesh;
const std::vector<IO::MeshDatabase> database = d_database[timestate];
const std::string timestep = d_timesteps[timestate];
for (size_t i=0; i<database.size(); i++) {
if ( database[i].name==std::string(meshname) ) {
DebugStream2 << " calling getMesh" << std::endl;
try {
mesh = IO::getMesh(d_path,timestep,database[i],domain);
} catch (const std::exception &err) {
DebugStream1 << " Caught errror calling getMesh:" << std::endl;
DebugStream1 << err.what() << std::endl;
} catch (...) {
DebugStream1 << " Caught unknown errror calling getMesh" << std::endl;
return NULL;
}
}
}
if ( mesh==NULL ) {
DebugStream1 << " Error loading mesh" << std::endl;
return NULL;
}
// Create the mesh in vtk
vtkDataSet* vtkMesh = meshToVTK(mesh);
DebugStream2 << " mesh created:" << std::endl;
ASSERT(vtkMesh!=NULL);
DebugStream2 << " " << vtkMesh->GetNumberOfCells() << std::endl;
vtkMesh->PrintSelf(DebugStream2,vtkIndent(6));
// Cache the mesh and return
// meshcache[cache_name] = mesh;
get_time(&stop);
DebugStream2 << " Time required: " << get_diff(start,stop,freq) << std::endl;
return vtkMesh;
}
// ****************************************************************************
// Method: avtLBMFileFormat::GetVar
//
// Purpose:
// Gets a scalar variable associated with this file. Although VTK has
// support for many different types, the best bet is vtkFloatArray, since
// that is supported everywhere through VisIt.
//
// Arguments:
// timestate The index of the timestate. If GetNTimesteps returned
// 'N' time steps, this is guaranteed to be between 0 and N-1.
// domain The index of the domain. If there are NDomains, this
// value is guaranteed to be between 0 and NDomains-1,
// regardless of block origin.
// varname The name of the variable requested.
//
// Programmer: mbt -- generated by xml2avt
// Creation: Wed Jul 9 10:59:08 PDT 2014
//
// ****************************************************************************
vtkDataArray *
avtLBMFileFormat::GetVar(int timestate, int domain, const char *meshvarname)
{
DebugStream1 << "avtLBMFileFormat::GetVar: " << meshvarname
<< "," << timestate << "," << domain << std::endl;
std::vector<std::string> tmp = IO::splitList(meshvarname,'/');
ASSERT(tmp.size()==2);
std::string varname = tmp[0];
std::string meshname = tmp[1];
const std::vector<IO::MeshDatabase> database = d_database[timestate];
const std::string timestep = d_timesteps[timestate];
std::shared_ptr<const IO::Variable> variable;
for (size_t i=0; i<database.size(); i++) {
if ( database[i].name==std::string(meshname) ) {
DebugStream2 << " calling getVar" << std::endl;
try {
variable = IO::getVariable(d_path,timestep,database[i],domain,varname);
} catch (const std::exception &err) {
DebugStream1 << " Caught errror calling getVar:" << std::endl;
DebugStream1 << err.what() << std::endl;
} catch (...) {
DebugStream1 << " Caught unknown errror calling getVar" << std::endl;
return NULL;
}
}
}
if ( variable == NULL )
EXCEPTION1(InvalidVariableException, varname);
vtkDataArray* vtkVar = varToVTK(variable);
return vtkVar;
}
// ****************************************************************************
// Method: avtLBMFileFormat::GetVectorVar
//
// Purpose:
// Gets a vector variable associated with this file. Although VTK has
// support for many different types, the best bet is vtkFloatArray, since
// that is supported everywhere through VisIt.
//
// Arguments:
// timestate The index of the timestate. If GetNTimesteps returned
// 'N' time steps, this is guaranteed to be between 0 and N-1.
// domain The index of the domain. If there are NDomains, this
// value is guaranteed to be between 0 and NDomains-1,
// regardless of block origin.
// varname The name of the variable requested.
//
// Programmer: mbt -- generated by xml2avt
// Creation: Wed Jul 9 10:59:08 PDT 2014
//
// ****************************************************************************
vtkDataArray *
avtLBMFileFormat::GetVectorVar(int timestate, int domain, const char *varname)
{
cerr << "avtLBMFileFormat::GetVectorVar - " << varname
<< "," << timestate << "," << domain << std::endl;
DebugStream1 << "avtLBMFileFormat::GetVectorVar" << std::endl;
EXCEPTION1(InvalidVariableException, varname);
return NULL;
//
// If you have a file format where variables don't apply (for example a
// strictly polygonal format like the STL (Stereo Lithography) format,
// then uncomment the code below.
//
// EXCEPTION1(InvalidVariableException, varname);
//
//
// If you do have a vector variable, here is some code that may be helpful.
//
// int ncomps = YYY; // This is the rank of the vector - typically 2 or 3.
// int ntuples = XXX; // this is the number of entries in the variable.
// vtkFloatArray *rv = vtkFloatArray::New();
// int ucomps = (ncomps == 2 ? 3 : ncomps);
// rv->SetNumberOfComponents(ucomps);
// rv->SetNumberOfTuples(ntuples);
// float *one_entry = new float[ucomps];
// for (int i = 0 ; i < ntuples ; i++)
// {
// int j;
// for (j = 0 ; j < ncomps ; j++)
// one_entry[j] = ...
// for (j = ncomps ; j < ucomps ; j++)
// one_entry[j] = 0.;
// rv->SetTuple(i, one_entry);
// }
//
// delete [] one_entry;
// return rv;
//
}