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

Adding test for blob identification. Adding visit capabilities for domain and volume variables

Browse Source
This commit is contained in:
Mark Berrill
2015-06-11 15:01:24 -04:00
parent db1650ea5d
commit 5fc08ce163
18 changed files with 558 additions and 60 deletions
Download Patch File Download Diff File Expand all files Collapse all files

66
IO/Mesh.cpp
View File

@@ -291,6 +291,72 @@ void TriMesh::unpack( const std::pair<size_t,void*>& data_in )
}
/****************************************************
* Domain mesh *
****************************************************/
DomainMesh::DomainMesh():
nprocx(0), nprocy(0), nprocz(0), rank(0),
nx(0), ny(0), nz(0),
Lx(0), Ly(0), Lz(0)
{
}
DomainMesh::DomainMesh( RankInfoStruct data,
int nx2, int ny2, int nz2, double Lx2, double Ly2, double Lz2 ):
nprocx(data.nx), nprocy(data.ny), nprocz(data.nz), rank(data.rank[1][1][1]),
nx(nx2), ny(ny2), nz(nz2),
Lx(Lx2), Ly(Ly2), Lz(Lz2)
{
}
DomainMesh::~DomainMesh()
{
}
size_t DomainMesh::numberPointsVar( VariableType type ) const
{
size_t N = 0;
if ( type==NodeVariable )
N = (nx+1)*(ny+1)*(nz+1);
else if ( type==SurfaceVariable )
N = (nx+1)*ny*nz + nx*(ny+1)*nz + nx*ny*(nz+1);
else if ( type==VolumeVariable )
N = nx*ny*nz;
return N;
}
std::pair<size_t,void*> DomainMesh::pack( int level ) const
{
std::pair<size_t,void*> data(0,NULL);
data.first = 7*sizeof(double);
data.second = new double[7];
int *data_int = reinterpret_cast<int*>(data.second);
double *data_double = &reinterpret_cast<double*>(data.second)[4];
data_int[0] = nprocx;
data_int[1] = nprocy;
data_int[2] = nprocz;
data_int[3] = rank;
data_int[4] = nx;
data_int[5] = ny;
data_int[6] = nz;
data_double[0] = Lx;
data_double[1] = Ly;
data_double[2] = Lz;
return data;
}
void DomainMesh::unpack( const std::pair<size_t,void*>& data )
{
const int *data_int = reinterpret_cast<const int*>(data.second);
const double *data_double = &reinterpret_cast<const double*>(data.second)[4];
nprocx = data_int[0];
nprocy = data_int[1];
nprocz = data_int[2];
rank = data_int[3];
nx = data_int[4];
ny = data_int[5];
nz = data_int[6];
Lx = data_double[0];
Ly = data_double[1];
Lz = data_double[2];
}
/****************************************************
* Converters *
****************************************************/

33
IO/Mesh.h
View File

@@ -5,7 +5,9 @@
#include <string.h>
#include <vector>
#include "common/Array.h"
#include "common/PointList.h"
#include "common/Communication.h"
#include "shared_ptr.h"
@@ -129,9 +131,36 @@ public:
};
/*! \class Domain
\brief A class used to hold the domain
*/
class DomainMesh: public Mesh
{
public:
//! Empty constructor
DomainMesh();
//! Default constructor
DomainMesh( RankInfoStruct rank_data, int nx, int ny, int nz, double Lx, double Ly, double Lz );
//! Destructor
virtual ~DomainMesh();
//! Mesh class name
virtual std::string className() const { return "DomainMesh"; }
//! Number of points for the given variable type
virtual size_t numberPointsVar( VariableType type ) const;
//! Pack the data
virtual std::pair<size_t,void*> pack( int level ) const;
//! Unpack the data
virtual void unpack( const std::pair<size_t,void*>& data );
public:
int nprocx, nprocy, nprocz, rank;
int nx, ny, nz;
double Lx, Ly, Lz;
};
/*! \class Variable
\brief A base class fore variables
\brief A base class for variables
*/
struct Variable
{
@@ -140,7 +169,7 @@ public:
unsigned int dim; //!< Number of points per grid point (1: scalar, 3: vector, ...)
VariableType type; //!< Variable type
std::string name; //!< Variable name
std::vector<double> data; //!< Variable data
Array<double> data; //!< Variable data
//! Empty constructor
Variable(): type(NullVariable) {}
//! Destructor

10
IO/MeshDatabase.cpp
View File

@@ -399,14 +399,16 @@ std::vector<MeshDatabase> read( const std::string& filename )
// Return the mesh type
IO::MeshType meshType( std::shared_ptr<IO::Mesh> mesh )
IO::MeshType meshType( const IO::Mesh& mesh )
{
IO::MeshType type = IO::Unknown;
if ( std::dynamic_pointer_cast<IO::PointList>(mesh).get()!=NULL ) {
const std::string meshClass = mesh.className();
if ( meshClass=="PointList" ) {
type = IO::PointMesh;
} else if ( std::dynamic_pointer_cast<IO::TriList>(mesh).get()!=NULL ||
std::dynamic_pointer_cast<IO::TriMesh>(mesh).get()!=NULL ) {
} else if ( meshClass=="TriList" || meshClass=="TriMesh" ) {
type = IO::SurfaceMesh;
} else if ( meshClass=="DomainMesh" ) {
type = IO::VolumeMesh;
} else {
ERROR("Unknown mesh");
}

2
IO/MeshDatabase.h
View File

@@ -81,7 +81,7 @@ std::vector<MeshDatabase> read( const std::string& filename );
//! Return the mesh type
IO::MeshType meshType( std::shared_ptr<IO::Mesh> mesh );
IO::MeshType meshType( const IO::Mesh& mesh );
} // IO namespace

10
IO/Reader.cpp
View File

@@ -100,6 +100,9 @@ std::shared_ptr<IO::Mesh> IO::getMesh( const std::string& path, const std::strin
C[i].z = data[9*i+8];
}
mesh = trilist;
} else if ( meshDatabase.type==IO::VolumeMesh ) {
// this was never supported in the old format
mesh = std::shared_ptr<DomainMesh>( new DomainMesh() );
} else {
ERROR("Unknown mesh type");
}
@@ -124,6 +127,8 @@ std::shared_ptr<IO::Mesh> IO::getMesh( const std::string& path, const std::strin
mesh.reset( new IO::TriMesh() );
} else if ( meshDatabase.meshClass=="TriList" ) {
mesh.reset( new IO::TriList() );
} else if ( meshDatabase.meshClass=="DomainMesh" ) {
mesh.reset( new IO::DomainMesh() );
} else {
ERROR("Unknown mesh class");
}
@@ -170,11 +175,8 @@ std::shared_ptr<IO::Variable> IO::getVariable( const std::string& path, const st
var->type = static_cast<IO::VariableType>(type);
var->name = variable;
var->data.resize(N);
double *var_data = NULL;
if ( !var->data.empty() )
var_data = &var->data[0];
if ( precision=="double" ) {
memcpy(var_data,data,bytes);
memcpy(var->data.get(),data,bytes);
} else {
ERROR("Format not implimented");
}

16
IO/Writer.cpp
View File

@@ -29,7 +29,7 @@ static std::vector<IO::MeshDatabase> writeMeshesOrigFormat( const std::vector<IO
std::shared_ptr<IO::Mesh> mesh = meshData[i].mesh;
IO::MeshDatabase mesh_entry;
mesh_entry.name = meshData[i].meshName;
mesh_entry.type = meshType(mesh);
mesh_entry.type = meshType(*mesh);
mesh_entry.meshClass = meshData[i].mesh->className();
mesh_entry.format = 1;
IO::DatabaseEntry domain;
@@ -42,7 +42,8 @@ static std::vector<IO::MeshDatabase> writeMeshesOrigFormat( const std::vector<IO
//for (size_t j=0; j<meshData[i].vars.size(); j++)
// mesh_entry.variables.push_back( meshData[i].vars[j]->name );
}
if ( std::dynamic_pointer_cast<IO::PointList>(mesh).get()!=NULL ) {
const std::string meshClass = mesh->className();
if ( meshClass=="PointList" ) {
// List of points
std::shared_ptr<IO::PointList> pointlist = std::dynamic_pointer_cast<IO::PointList>(mesh);
const std::vector<Point>& P = pointlist->points;
@@ -51,8 +52,7 @@ static std::vector<IO::MeshDatabase> writeMeshesOrigFormat( const std::vector<IO
x[0] = P[i].x; x[1] = P[i].y; x[2] = P[i].z;
fwrite(x,sizeof(double),3,fid);
}
} else if ( std::dynamic_pointer_cast<IO::TriList>(mesh).get()!=NULL ||
std::dynamic_pointer_cast<IO::TriMesh>(mesh).get()!=NULL ) {
} else if ( meshClass=="TriList" || meshClass=="TriMesh" ) {
// Triangle mesh
std::shared_ptr<IO::TriList> trilist = IO::getTriList(mesh);
const std::vector<Point>& A = trilist->A;
@@ -65,6 +65,8 @@ static std::vector<IO::MeshDatabase> writeMeshesOrigFormat( const std::vector<IO
tri[6] = C[i].x; tri[7] = C[i].y; tri[8] = C[i].z;
fwrite(tri,sizeof(double),9,fid);
}
} else if ( meshClass=="DomainMesh" ) {
// This format was never supported with the old format
} else {
ERROR("Unknown mesh");
}
@@ -87,7 +89,7 @@ static IO::MeshDatabase write_domain( FILE *fid, const std::string& filename,
// Create the MeshDatabase
IO::MeshDatabase database;
database.name = mesh.meshName;
database.type = meshType(mesh.mesh);
database.type = meshType(*(mesh.mesh));
database.meshClass = mesh.mesh->className();
database.format = format;
// Write the mesh
@@ -117,8 +119,8 @@ static IO::MeshDatabase write_domain( FILE *fid, const std::string& filename,
std::pair<std::pair<std::string,std::string>,IO::DatabaseEntry>(key,variable) );
int dim = mesh.vars[i]->dim;
int type = static_cast<int>(mesh.vars[i]->type);
size_t N = mesh.vars[i]->data.size();
const void* data = N==0 ? 0:&mesh.vars[i]->data[0];
size_t N = mesh.vars[i]->data.length();
const double* data = N==0 ? NULL:mesh.vars[i]->data.get();
if ( type == static_cast<int>(IO::NullVariable) ) {
ERROR("Variable type not set");
}

2
cmake/Find_TIMER.cmake
View File

@@ -32,7 +32,7 @@ FUNCTION( CONFIGURE_TIMER DEFAULT_USE_TIMER NULL_TIMER_DIR )
VERIFY_PATH( ${TIMER_DIRECTORY}/lib )
FIND_LIBRARY( TIMER_LIBS NAMES timerutility PATHS ${TIMER_DIRECTORY}/lib NO_DEFAULT_PATH )
SET( TIMER_INCLUDE ${TIMER_DIRECTORY}/include )
SET( TIMER_CXXFLAGS "-DUSE_TIMER -I${TIMER_DIRECTORY}/include" )
SET( TIMER_CXXFLAGS -DUSE_TIMER )
SET( TIMER_LDFLAGS -L${TIMER_DIRECTORY}/lib )
SET( TIMER_LDLIBS -ltimerutility )
ELSE()

10
common/Communication.h
View File

@@ -60,11 +60,19 @@ public:
/*!
* @brief Communicate the halos
* @param[in] array The array on which we fill the halos
*/
void fill( Array<TYPE>& array );
/*!
* @brief Copy data from the src array to the dst array
* @param[in] src The src array with or without halos
* @param[in] dst The dst array with or without halos
*/
template<class TYPE1, class TYPE2>
void copy( const Array<TYPE1>& src, Array<TYPE2>& dst );
private:
RankInfoStruct info;
int nx, ny, nz, ngx, ngy, ngz, depth;

56
common/Communication.hpp
View File

@@ -199,4 +199,60 @@ void fillHalo<TYPE>::unpack( Array<TYPE>& data, int i0, int j0, int k0, const TY
}
/********************************************************
* Function to remove the ghost halo *
********************************************************/
template<class TYPE>
template<class TYPE1, class TYPE2>
void fillHalo<TYPE>::copy( const Array<TYPE1>& src, Array<TYPE2>& dst )
{
ASSERT(src.size(0)==nx||src.size(0)==nx+2*ngx);
ASSERT(dst.size(0)==nx||dst.size(0)==nx+2*ngx);
bool src_halo = src.size(0)==nx+2*ngx;
bool dst_halo = dst.size(0)==nx+2*ngx;
if ( src_halo ) {
ASSERT(src.size(0)==nx+2*ngx);
ASSERT(src.size(1)==ny+2*ngy);
ASSERT(src.size(2)==nz+2*ngz);
} else {
ASSERT(src.size(0)==nx);
ASSERT(src.size(1)==ny);
ASSERT(src.size(2)==nz);
}
if ( dst_halo ) {
ASSERT(dst.size(0)==nx+2*ngx);
ASSERT(dst.size(1)==ny+2*ngy);
ASSERT(dst.size(2)==nz+2*ngz);
} else {
ASSERT(dst.size(0)==nx);
ASSERT(dst.size(1)==ny);
ASSERT(dst.size(2)==nz);
}
if ( src_halo == dst_halo ) {
// Src and dst halos match
for (size_t i=0; i<src.length(); i++)
dst(i) = src(i);
} else if ( src_halo && !dst_halo ) {
// Src has halos
for (size_t k=0; k<nz; k++) {
for (size_t j=0; j<ny; j++) {
for (size_t i=0; i<nx; i++) {
dst(i,j,k) = src(i+ngx,j+ngy,k+ngz);
}
}
}
} else if ( !src_halo && dst_halo ) {
// Src has halos
for (size_t k=0; k<nz; k++) {
for (size_t j=0; j<ny; j++) {
for (size_t i=0; i<nx; i++) {
dst(i+ngx,j+ngy,k+ngz) = src(i,j,k);
}
}
}
fill(dst);
}
}
#endif

51
common/Utilities.cpp
View File

@@ -7,6 +7,8 @@
#include <fstream>
#include <string.h>
#include <signal.h>
#include <math.h>
#include <algorithm>
#ifdef USE_MPI
#include "mpi.h"
@@ -388,3 +390,52 @@ std::vector<std::string> Utilities::getCallStack()
#error Unknown OS
#endif
// Factor a number into it's prime factors
std::vector<int> Utilities::factor(size_t number)
{
if ( number<=3 )
return std::vector<int>(1,(int)number);
size_t i, n, n_max;
bool factor_found;
// Compute the maximum number of factors
int N_primes_max = 1;
n = number;
while (n >>= 1) ++N_primes_max;
// Initialize n, factors
n = number;
std::vector<int> factors;
factors.reserve(N_primes_max);
while ( 1 ) {
// Check if n is a trivial prime number
if ( n==2 || n==3 || n==5 ) {
factors.push_back( (int) n );
break;
}
// Check if n is divisible by 2
if ( n%2 == 0 ) {
factors.push_back( 2 );
n/=2;
continue;
}
// Check each odd number until a factor is reached
n_max = (size_t) floor(sqrt((double) n));
factor_found = false;
for (i=3; i<=n_max; i+=2) {
if ( n%i == 0 ) {
factors.push_back( i );
n/=i;
factor_found = true;
break;
}
}
if ( factor_found )
continue;
// No factors were found, the number must be prime
factors.push_back( (int) n );
break;
}
// Sort the factors
std::sort( factors.begin(), factors.end() );
return factors;
}

3
common/Utilities.h
View File

@@ -60,6 +60,9 @@ namespace Utilities
//! Function to get the resolution of time
double tick();
//! Factor a number into it's prime factors
std::vector<int> factor(size_t number);
} // namespace Utilities

1
tests/CMakeLists.txt
View File

@@ -25,6 +25,7 @@ ADD_LBPM_TEST( TestContactAngle )
ADD_LBPM_TEST_1_2_4( TestTwoPhase )
ADD_LBPM_TEST_PARALLEL( TestTwoPhase 8 )
ADD_LBPM_TEST_PARALLEL( TestBlobAnalyze 8 )
ADD_LBPM_TEST_PARALLEL( TestBlobIdentify 1 )
ADD_LBPM_TEST_PARALLEL( TestSegDist 8 )
ADD_LBPM_TEST_1_2_4( testCommunication )
ADD_LBPM_TEST_1_2_4( testUtilities )

185
tests/TestBlobIdentify.cpp Normal file
View File

@@ -0,0 +1,185 @@
// Sequential blob analysis
// Reads parallel simulation data and performs connectivity analysis
// and averaging on a blob-by-blob basis
// James E. McClure 2014
#include <iostream>
#include <math.h>
#include "common/pmmc.h"
#include "common/Communication.h"
#include "analysis/analysis.h"
#include "IO/MeshDatabase.h"
#include "IO/Reader.h"
#include "IO/Writer.h"
#include "ProfilerApp.h"
// Get a random number in [0 1]
inline double rand2()
{
return static_cast<double>(rand())/static_cast<double>(RAND_MAX);
}
struct bubble_struct {
Point center;
double radius;
// Get the distance to the bubble
double dist( const Point& p, double Lx, double Ly, double Lz ) {
double x = std::min(fabs(p.x-center.x),std::min(fabs(p.x-center.x-Lx),fabs(p.x-center.x+Lx)));
double y = std::min(fabs(p.y-center.y),std::min(fabs(p.y-center.y-Ly),fabs(p.y-center.y+Ly)));
double z = std::min(fabs(p.z-center.z),std::min(fabs(p.z-center.z-Lz),fabs(p.z-center.z+Lz)));
return sqrt(x*x+y*y+z*z)-radius;
}
// Check if this bubble overlaps with rhs
bool overlap( const bubble_struct& rhs, double Lx, double Ly, double Lz ) {
return dist(rhs.center,Lx,Ly,Lz) <= radius+rhs.radius;
}
// Create a random bubble
static bubble_struct rand( double Lx, double Ly, double Lz, double R0 ) {
bubble_struct bubble;
bubble.center.x = Lx*rand2();
bubble.center.y = Ly*rand2();
bubble.center.z = Lz*rand2();
bubble.radius = R0;
bubble.radius *= 1 + 0.4*(1.0-2.0*rand2());
return bubble;
}
};
// Create a random set of bubles
std::vector<bubble_struct> create_bubbles( int N_bubbles, double Lx, double Ly, double Lz )
{
int rank;
MPI_Comm_rank(MPI_COMM_WORLD,&rank);
std::vector<bubble_struct> bubbles(N_bubbles);
if ( rank == 0 ) {
double R0 = 0.2*Lx*Ly*Lz/pow((double)N_bubbles,0.333);
for (int i=0; i<N_bubbles; i++) {
bool overlap = true;
while ( overlap ) {
bubbles[i] = bubble_struct::rand(Lx,Ly,Lz,R0);
overlap = false;
for (int k=0; k<i; k++)
overlap = overlap || bubbles[i].overlap(bubbles[k],Lx,Ly,Lz);
}
}
}
size_t N_bytes = N_bubbles*sizeof(bubble_struct);
MPI_Bcast(&bubbles[0],N_bytes,MPI_CHAR,0,MPI_COMM_WORLD);
return bubbles;
}
// Main
int main(int argc, char **argv)
{
// Initialize MPI
int rank, nprocs;
MPI_Init(&argc,&argv);
MPI_Comm_rank(MPI_COMM_WORLD,&rank);
MPI_Comm_size(MPI_COMM_WORLD,&nprocs);
PROFILE_ENABLE(0);
PROFILE_DISABLE_TRACE();
PROFILE_SYNCHRONIZE();
PROFILE_START("main");
if ( rank==0 ) {
printf("-----------------------------------------------------------\n");
printf("Testing Blob Identification \n");
printf("-----------------------------------------------------------\n");
}
// Set the domain information
std::vector<int> factors = Utilities::factor(nprocs);
int nproc[3]={1,1,1};
for (size_t i=0; i<factors.size(); i++)
nproc[i%3] *= factors[i];
int nx = 100/nproc[0];
int ny = 120/nproc[1];
int nz = 110/nproc[2];
double Lx=1.0, Ly=1.0, Lz=1.0;
// Get the rank info
const RankInfoStruct rank_info(rank,nproc[0],nproc[1],nproc[2]);
// Create the dummy info
DoubleArray Phase(nx+2,ny+2,nz+2);
DoubleArray SignDist(nx+2,ny+2,nz+2);
Phase.fill(1);
SignDist(0);
std::vector<bubble_struct> bubbles = create_bubbles(20,Lx,Ly,Lz);
for (int k=0; k<nz; k++) {
double z = Lz*(nz*rank_info.kz+k+0.5)/(nz*nproc[2]);
for (int j=0; j<ny; j++) {
double y = Ly*(ny*rank_info.jy+j+0.5)/(ny*nproc[1]);
for (int i=0; i<nx; i++) {
double x = Lx*(nx*rank_info.ix+i+0.5)/(nx*nproc[0]);
double dist = 1e100;
for (size_t b=0; b<bubbles.size(); b++)
dist = std::min(dist,bubbles[b].dist(Point(x,y,z),Lx,Ly,Lz));
SignDist(i+1,j+1,k+1) = -dist;
}
}
}
// Communication the halos
fillHalo<double> fillData(rank_info,nx,ny,nz,1,1,1,0,1);
fillData.fill(Phase);
fillData.fill(SignDist);
// Find blob domains
if ( rank==0 ) { printf("Finding blob domains\n"); }
double vF=0.0;
double vS=0.0;
IntArray GlobalBlobID;
int nblobs = ComputeGlobalBlobIDs(nx,ny,nz,rank_info,
Phase,SignDist,vF,vS,GlobalBlobID);
if ( rank==0 ) { printf("Identified %i blobs\n",nblobs); }
// Create the MeshDataStruct
std::vector<IO::MeshDataStruct> meshData(1);
meshData[0].meshName = "domain";
meshData[0].mesh = std::shared_ptr<IO::DomainMesh>( new IO::DomainMesh(rank_info,nx,ny,nz,Lx,Ly,Lz) );
std::shared_ptr<IO::Variable> PhaseVar( new IO::Variable() );
std::shared_ptr<IO::Variable> SignDistVar( new IO::Variable() );
std::shared_ptr<IO::Variable> BlobIDVar( new IO::Variable() );
PhaseVar->name = "phase";
PhaseVar->type = IO::VolumeVariable;
PhaseVar->dim = 1;
PhaseVar->data.resize(nx,ny,nz);
meshData[0].vars.push_back(PhaseVar);
SignDistVar->name = "SignDist";
SignDistVar->type = IO::VolumeVariable;
SignDistVar->dim = 1;
SignDistVar->data.resize(nx,ny,nz);
meshData[0].vars.push_back(SignDistVar);
BlobIDVar->name = "BlobID";
BlobIDVar->type = IO::VolumeVariable;
BlobIDVar->dim = 1;
BlobIDVar->data.resize(nx,ny,nz);
meshData[0].vars.push_back(BlobIDVar);
// Save the results
fillData.copy(Phase,PhaseVar->data);
fillData.copy(SignDist,SignDistVar->data);
fillData.copy(GlobalBlobID,BlobIDVar->data);
IO::writeData( 0, meshData, 2 );
// Check the results
int N_errors = 0;
if ( nblobs != (int) bubbles.size() ) {
printf("Error, detected number of bubbles %i does not match expected %i\n",nblobs,(int)bubbles.size());
N_errors++;
}
// Finished
PROFILE_STOP("main");
PROFILE_SAVE("TestBlobIdentify",false);
MPI_Barrier(MPI_COMM_WORLD);
MPI_Finalize();
return 0;
}

8
tests/TestBubble.cpp
View File

@@ -2170,14 +2170,14 @@ int main(int argc, char **argv)
dist->name = "distance";
dist->dim = 1;
dist->type = IO::NodeVariable;
dist->data.resize(3*mesh->A.size());
dist->data.resize(3,mesh->A.size());
for (size_t i=0; i<mesh->A.size(); i++) {
const Point& a = mesh->A[i];
const Point& b = mesh->B[i];
const Point& c = mesh->C[i];
dist->data[3*i+0] = sqrt(a.x*a.x+a.y*a.y+a.z*a.z);
dist->data[3*i+1] = sqrt(b.x*b.x+b.y*b.y+b.z*b.z);
dist->data[3*i+2] = sqrt(c.x*c.x+c.y*c.y+c.z*c.z);
dist->data(0,i) = sqrt(a.x*a.x+a.y*a.y+a.z*a.z);
dist->data(1,i) = sqrt(b.x*b.x+b.y*b.y+b.z*b.z);
dist->data(2,i) = sqrt(c.x*c.x+c.y*c.y+c.z*c.z);
}
meshData[k].vars.push_back(dist);
}

66
tests/TestWriter.cpp
View File

@@ -78,28 +78,42 @@ int main(int argc, char **argv)
return -1;
}
}
RankInfoStruct rank_data( rank, nprocs, 1, 1 );
std::shared_ptr<IO::DomainMesh> domain( new IO::DomainMesh(rank_data,6,7,8,1.0,1.0,1.0) );
// Create the variables
std::shared_ptr<IO::Variable> dist_set1( new IO::Variable() );
std::shared_ptr<IO::Variable> dist_list( new IO::Variable() );
std::shared_ptr<IO::Variable> domain_var( new IO::Variable() );
dist_set1->dim = 1;
dist_list->dim = 1;
domain_var->dim = 1;
dist_set1->name = "Distance";
dist_list->name = "Distance";
domain_var->name = "Distance";
dist_set1->type = IO::NodeVariable;
dist_list->type = IO::NodeVariable;
domain_var->type = IO::VolumeVariable;
dist_set1->data.resize( N_points );
for (int i=0; i<N_points; i++)
dist_set1->data[i] = distance(set1->points[i]);
dist_list->data.resize( 3*N_tri );
dist_set1->data(i) = distance(set1->points[i]);
dist_list->data.resize( 3, N_tri );
for (int i=0; i<N_tri; i++) {
dist_list->data[3*i+0] = distance(trilist->A[i]);
dist_list->data[3*i+1] = distance(trilist->B[i]);
dist_list->data[3*i+2] = distance(trilist->C[i]);
dist_list->data(0,i) = distance(trilist->A[i]);
dist_list->data(1,i) = distance(trilist->B[i]);
dist_list->data(2,i) = distance(trilist->C[i]);
}
domain_var->data.resize(domain->nx,domain->ny,domain->nz);
for (int i=0; i<domain->nx; i++) {
for (int j=0; j<domain->ny; j++) {
for (int k=0; k<domain->nz; k++) {
domain_var->data(i,j,k) = distance(Point(i,j,k));
}
}
}
// Create the MeshDataStruct
std::vector<IO::MeshDataStruct> meshData(3);
std::vector<IO::MeshDataStruct> meshData(4);
meshData[0].meshName = "pointmesh";
meshData[0].mesh = set1;
meshData[0].vars.push_back(dist_set1);
@@ -109,10 +123,16 @@ int main(int argc, char **argv)
meshData[2].meshName = "trilist";
meshData[2].mesh = trilist;
meshData[2].vars.push_back(dist_list);
meshData[3].meshName = "domain";
meshData[3].mesh = domain;
meshData[3].vars.push_back(domain_var);
// Write the data
IO::writeData( 0, meshData, 1 );
IO::writeData( 3, meshData, 2 );
std::vector<int> format(2,0);
format[0] = 2;
format[1] = 1;
IO::writeData( 0, meshData, format[0] );
IO::writeData( 3, meshData, format[1] );
MPI_Barrier(MPI_COMM_WORLD);
// Get a list of the timesteps
@@ -142,14 +162,14 @@ int main(int argc, char **argv)
}
// For each domain, load the mesh and check its data
for (size_t j=0; j<list.size(); j++) {
for (size_t k=0; k<list[i].domains.size(); k++) {
for (size_t k=0; k<list[j].domains.size(); k++) {
std::shared_ptr<IO::Mesh> mesh = IO::getMesh(".",timesteps[i],list[j],k);
if ( mesh.get()==NULL ) {
printf("Failed to load %s\n",meshData[i].meshName.c_str());
printf("Failed to load %s\n",list[j].name.c_str());
pass = false;
break;
}
if ( meshData[j].meshName=="pointmesh" ) {
if ( list[j].name=="pointmesh" ) {
// Check the pointmesh
std::shared_ptr<IO::PointList> pmesh = IO::getPointList(mesh);
if ( pmesh.get()==NULL ) {
@@ -161,7 +181,7 @@ int main(int argc, char **argv)
break;
}
}
if ( meshData[j].meshName=="trimesh" || meshData[j].meshName=="trilist" ) {
if ( list[j].name=="trimesh" || list[j].name=="trilist" ) {
// Check the trimesh/trilist
std::shared_ptr<IO::TriMesh> mesh1 = IO::getTriMesh(mesh);
std::shared_ptr<IO::TriList> mesh2 = IO::getTriList(mesh);
@@ -192,6 +212,16 @@ int main(int argc, char **argv)
pass = false;
}
}
if ( list[j].name=="domain" && format[i]>2 ) {
// Check the domain mesh
const IO::DomainMesh& mesh1 = *std::dynamic_pointer_cast<IO::DomainMesh>(mesh);
if ( mesh1.nprocx!=domain->nprocx || mesh1.nprocy!=domain->nprocy || mesh1.nprocz!=domain->nprocz )
pass = false;
if ( mesh1.nx!=domain->nx || mesh1.ny!=domain->ny || mesh1.nz!=domain->nz )
pass = false;
if ( mesh1.Lx!=domain->Lx || mesh1.Ly!=domain->Ly || mesh1.Lz!=domain->Lz )
pass = false;
}
}
}
if ( pass ) {
@@ -201,7 +231,7 @@ int main(int argc, char **argv)
continue;
}
// For each domain, load the variables and check their data
if ( i==0 )
if ( format[i]==1 )
continue; // Format 1 does not support variables
for (size_t j=0; j<list.size(); j++) {
const IO::MeshDataStruct* mesh0 = NULL;
@@ -211,8 +241,8 @@ int main(int argc, char **argv)
break;
}
}
for (size_t v=0; v<list[i].variables.size(); v++) {
for (size_t k=0; k<list[i].domains.size(); k++) {
for (size_t v=0; v<list[j].variables.size(); v++) {
for (size_t k=0; k<list[j].domains.size(); k++) {
std::shared_ptr<const IO::Variable> variable =
IO::getVariable(".",timesteps[i],list[j],k,list[j].variables[v].name);
const IO::Variable& var1 = *mesh0->vars[v];
@@ -220,10 +250,10 @@ int main(int argc, char **argv)
pass = var1.name == var2.name;
pass = pass && var1.dim == var2.dim;
pass = pass && var1.type == var2.type;
pass = pass && var1.data.size() == var2.data.size();
pass = pass && var1.data.length() == var2.data.length();
if ( pass ) {
for (size_t m=0; m<var1.data.size(); m++)
pass = pass && approx_equal(var1.data[m],var2.data[m]);
for (size_t m=0; m<var1.data.length(); m++)
pass = pass && approx_equal(var1.data(m),var2.data(m));
}
if ( pass ) {
ut.passes("Variables match");

1
visit/LBM.xml
View File

@@ -5,6 +5,7 @@
</FilePatterns>
<CXXFLAGS>
-I${LBPM_INSTALL_DIR}/include
-I${TIMER_INCLUDE}
${MPI_CXXFLAGS}
${TIMER_CXXFLAGS}
</CXXFLAGS>

8
visit/avtLBMFileFormat.C
View File

@@ -216,7 +216,7 @@ avtLBMFileFormat::PopulateDatabaseMetaData(avtDatabaseMetaData *md, int timestat
mmd->name = database[i].name;
mmd->meshType = vtkMeshType(database[i].type);
mmd->spatialDimension = 3;
mmd->topologicalDimension = -1;
mmd->topologicalDimension = vtkTopDim(database[i].type);
if ( mmd->meshType==AVT_SURFACE_MESH )
mmd->topologicalDimension = 2;
mmd->numBlocks = database[i].domains.size();
@@ -240,11 +240,11 @@ avtLBMFileFormat::PopulateDatabaseMetaData(avtDatabaseMetaData *md, int timestat
IO::VariableDatabase variable = database[i].variables[j];
std::string varname = variable.name + "/" + mmd->name;
avtCentering center = AVT_UNKNOWN_CENT;
if ( variable.type==IO::VariableType::NodeVariable ) {
if ( variable.type==IO::NodeVariable ) {
center = AVT_NODECENT;
} else if ( variable.type==IO::VariableType::SurfaceVariable ) {
} else if ( variable.type==IO::SurfaceVariable ) {
center = AVT_ZONECENT;
} else if ( variable.type==IO::VariableType::VolumeVariable ) {
} else if ( variable.type==IO::VolumeVariable ) {
center = AVT_ZONECENT;
}
if ( variable.dim==1 ) {

90
visit/vtkHelpers.h
View File

@@ -7,10 +7,17 @@
#include <vtkCellType.h>
#include <vtkPolyData.h>
#include <vtkCellArray.h>
#include <vtkDoubleArray.h>
#include <vtkStructuredGrid.h>
// LBPM headers
#include "IO/Mesh.h"
#include "IO/MeshDatabase.h"
extern std::ostream DebugStream1;
extern std::ostream DebugStream2;
// Convert a point array to vtkFloatArray
@@ -29,16 +36,30 @@ inline vtkFloatArray* pointToVTK( const std::vector<Point>& points )
inline avtMeshType vtkMeshType( IO::MeshType meshType )
{
avtMeshType vtkType = AVT_UNKNOWN_MESH;
if ( meshType==IO::MeshType::PointMesh )
if ( meshType==IO::PointMesh )
vtkType = AVT_POINT_MESH;
else if ( meshType==IO::MeshType::SurfaceMesh )
else if ( meshType==IO::SurfaceMesh )
vtkType = AVT_SURFACE_MESH;
else if ( meshType==IO::MeshType::VolumeMesh )
vtkType = AVT_UNSTRUCTURED_MESH;
else if ( meshType==IO::VolumeMesh )
vtkType = AVT_RECTILINEAR_MESH;
return vtkType;
}
// Return the topological dimension
inline int vtkTopDim( IO::MeshType meshType )
{
int dim = -1;
if ( meshType==IO::PointMesh )
dim = 1;
else if ( meshType==IO::SurfaceMesh )
dim = 2;
else if ( meshType==IO::VolumeMesh )
dim = 3;
return dim;
}
// Convert a PointList to vtkDataSet
inline vtkDataSet* PointListToVTK( std::shared_ptr<const IO::PointList> mesh )
{
@@ -119,22 +140,63 @@ inline vtkDataSet* TriListToVTK( std::shared_ptr<const IO::TriList> mesh )
}
// Convert a DomainMesh to vtkDataSet
inline vtkDataSet* DomainToVTK( std::shared_ptr<const IO::DomainMesh> mesh )
{
int nx = mesh->nx;
int ny = mesh->ny;
int nz = mesh->nz;
if ( nx==0 || ny==0 || nz==0 ) {
DebugStream::Stream1() << " Domain mesh is empty" << std::endl;
return NULL;
}
RankInfoStruct rank_data(mesh->rank,mesh->nprocx,mesh->nprocy,mesh->nprocz);
vtkFloatArray *x = vtkFloatArray::New();
vtkFloatArray *y = vtkFloatArray::New();
vtkFloatArray *z = vtkFloatArray::New();
for(int i=0; i<=nx; i++)
x->InsertNextValue(static_cast<double>(nx*rank_data.ix+i)/static_cast<double>(nx*mesh->nprocx));
for(int j=0; j<=ny; j++)
y->InsertNextValue(static_cast<double>(ny*rank_data.jy+j)/static_cast<double>(ny*mesh->nprocy));
for(int k=0; k<=nz; k++)
z->InsertNextValue(static_cast<double>(nz*rank_data.kz+k)/static_cast<double>(nz*mesh->nprocz));
vtkRectilinearGrid *vtkMesh = vtkRectilinearGrid::New();
vtkMesh->SetDimensions(nx+1,ny+1,nz+1);
vtkMesh->SetXCoordinates(x);
vtkMesh->SetYCoordinates(y);
vtkMesh->SetZCoordinates(z);
vtkMesh->ComputeBounds();
x->Delete();
y->Delete();
z->Delete();
return vtkMesh;
}
// Convert a mesh to vtkDataSet
inline vtkDataSet* meshToVTK( std::shared_ptr<const IO::Mesh> mesh )
{
if ( mesh.get() == NULL ){
DebugStream::Stream1() << " Mesh is NULL" << std::endl;
return NULL;
}
vtkDataSet* mesh2 = NULL;
if ( std::dynamic_pointer_cast<const IO::PointList>(mesh) != NULL ) {
const std::string className = mesh->className();
if ( className == "PointList" ) {
// We are dealing with a point mesh
mesh2 = PointListToVTK( std::dynamic_pointer_cast<const IO::PointList>(mesh) );
DebugStream::Stream1() << " Point mesh created" << std::endl;
} else if ( std::dynamic_pointer_cast<const IO::TriMesh>(mesh) != NULL ) {
DebugStream1 << " Point mesh created" << std::endl;
} else if ( className == "TriMesh" ) {
mesh2 = TriMeshToVTK( std::dynamic_pointer_cast<const IO::TriMesh>(mesh) );
DebugStream::Stream1() << " Surface mesh created" << std::endl;
} else if ( std::dynamic_pointer_cast<const IO::TriList>(mesh) != NULL ) {
DebugStream1 << " Surface mesh created" << std::endl;
} else if ( className == "TriList" ) {
mesh2 = TriListToVTK( std::dynamic_pointer_cast<const IO::TriList>(mesh) );
DebugStream::Stream1() << " Surface mesh created" << std::endl;
DebugStream1 << " Surface mesh created" << std::endl;
} else if ( className == "DomainMesh" ) {
mesh2 = DomainToVTK( std::dynamic_pointer_cast<const IO::DomainMesh>(mesh) );
DebugStream1 << " Volume mesh created" << std::endl;
} else {
DebugStream::Stream1() << " Error, unknown mesh type" << std::endl;
DebugStream1 << " Error, unknown mesh type" << std::endl;
return NULL;
}
return mesh2;
@@ -148,9 +210,9 @@ inline vtkDataArray* varToVTK( std::shared_ptr<const IO::Variable> var )
if ( var->dim==1 ) {
// Scalar variable
var2 = vtkFloatArray::New();
var2->SetNumberOfTuples(var->data.size());
for (size_t i=0; i<var->data.size(); i++)
var2->SetTuple1(i,var->data[i]);
var2->SetNumberOfTuples(var->data.length());
for (size_t i=0; i<var->data.length(); i++)
var2->SetTuple1(i,var->data(i));
} else {
DebugStream::Stream1() << " Error, variable not yet supported" << std::endl;
return NULL;