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

Replacing Array class with another Array class

Browse Source
This commit is contained in:
Mark Berrill 2015-04-22 14:23:55 -04:00
parent d33c1abd71
commit ef25839769
28 changed files with 1099 additions and 1050 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
CMakeLists.txt
View File

@ -104,7 +104,7 @@ IF ( NOT ONLY_BUILD_DOCS )
CONFIGURE_TIMER( 0 "${${PROJ}_INSTALL_DIR}/null_timer" )
CONFIGURE_LINE_COVERAGE()
INCLUDE( "${CMAKE_CURRENT_SOURCE_DIR}/cmake/SharedPtr.cmake" )
CONFIGURE_SHARED_PTR( "${LBPM_INSTALL_DIR}/include" "" )
CONFIGURE_SHARED_PTR( "${LBPM_INSTALL_DIR}/include" "std" )
ENDIF()

43
IO/Mesh.cpp
View File

@ -1,5 +1,6 @@
#include "Mesh.h"
#include "common/Utilities.h"
#include "shared_ptr.h"
#include <limits>
#include <stdint.h>
@ -192,7 +193,7 @@ TriMesh::TriMesh( size_t N_tri, size_t N_point )
B.resize(N_tri,-1);
C.resize(N_tri,-1);
}
TriMesh::TriMesh( size_t N_tri, shared_ptr<PointList> points )
TriMesh::TriMesh( size_t N_tri, std::shared_ptr<PointList> points )
{
vertices = points;
A.resize(N_tri,-1);
@ -293,45 +294,45 @@ void TriMesh::unpack( const std::pair<size_t,void*>& data_in )
/****************************************************
* Converters *
****************************************************/
shared_ptr<PointList> getPointList( shared_ptr<Mesh> mesh )
std::shared_ptr<PointList> getPointList( std::shared_ptr<Mesh> mesh )
{
return dynamic_pointer_cast<PointList>(mesh);
return std::dynamic_pointer_cast<PointList>(mesh);
}
shared_ptr<TriMesh> getTriMesh( shared_ptr<Mesh> mesh )
std::shared_ptr<TriMesh> getTriMesh( std::shared_ptr<Mesh> mesh )
{
shared_ptr<TriMesh> mesh2;
if ( dynamic_pointer_cast<TriMesh>(mesh).get() != NULL ) {
mesh2 = dynamic_pointer_cast<TriMesh>(mesh);
} else if ( dynamic_pointer_cast<TriList>(mesh).get() != NULL ) {
shared_ptr<TriList> trilist = dynamic_pointer_cast<TriList>(mesh);
std::shared_ptr<TriMesh> mesh2;
if ( std::dynamic_pointer_cast<TriMesh>(mesh).get() != NULL ) {
mesh2 = std::dynamic_pointer_cast<TriMesh>(mesh);
} else if ( std::dynamic_pointer_cast<TriList>(mesh).get() != NULL ) {
std::shared_ptr<TriList> trilist = std::dynamic_pointer_cast<TriList>(mesh);
ASSERT(trilist.get()!=NULL);
mesh2.reset( new TriMesh(*trilist) );
}
return mesh2;
}
shared_ptr<TriList> getTriList( shared_ptr<Mesh> mesh )
std::shared_ptr<TriList> getTriList( std::shared_ptr<Mesh> mesh )
{
shared_ptr<TriList> mesh2;
if ( dynamic_pointer_cast<TriList>(mesh).get() != NULL ) {
mesh2 = dynamic_pointer_cast<TriList>(mesh);
} else if ( dynamic_pointer_cast<TriMesh>(mesh).get() != NULL ) {
shared_ptr<TriMesh> trimesh = dynamic_pointer_cast<TriMesh>(mesh);
std::shared_ptr<TriList> mesh2;
if ( std::dynamic_pointer_cast<TriList>(mesh).get() != NULL ) {
mesh2 = std::dynamic_pointer_cast<TriList>(mesh);
} else if ( std::dynamic_pointer_cast<TriMesh>(mesh).get() != NULL ) {
std::shared_ptr<TriMesh> trimesh = std::dynamic_pointer_cast<TriMesh>(mesh);
ASSERT(trimesh.get()!=NULL);
mesh2.reset( new TriList(*trimesh) );
}
return mesh2;
}
shared_ptr<const PointList> getPointList( shared_ptr<const Mesh> mesh )
std::shared_ptr<const PointList> getPointList( std::shared_ptr<const Mesh> mesh )
{
return getPointList( const_pointer_cast<Mesh>(mesh) );
return getPointList( std::const_pointer_cast<Mesh>(mesh) );
}
shared_ptr<const TriMesh> getTriMesh( shared_ptr<const Mesh> mesh )
std::shared_ptr<const TriMesh> getTriMesh( std::shared_ptr<const Mesh> mesh )
{
return getTriMesh( const_pointer_cast<Mesh>(mesh) );
return getTriMesh( std::const_pointer_cast<Mesh>(mesh) );
}
shared_ptr<const TriList> getTriList( shared_ptr<const Mesh> mesh )
std::shared_ptr<const TriList> getTriList( std::shared_ptr<const Mesh> mesh )
{
return getTriList( const_pointer_cast<Mesh>(mesh) );
return getTriList( std::const_pointer_cast<Mesh>(mesh) );
}

20
IO/Mesh.h
View File

@ -108,7 +108,7 @@ public:
//! Constructor for Nt triangles and Np points
TriMesh( size_t N_tri, size_t N_point );
//! Constructor for Nt triangles and the given points
TriMesh( size_t N_tri, shared_ptr<PointList> points );
TriMesh( size_t N_tri, std::shared_ptr<PointList> points );
//! Constructor from TriList
TriMesh( const TriList& );
//! Destructor
@ -122,7 +122,7 @@ public:
//! Unpack the data
virtual void unpack( const std::pair<size_t,void*>& data );
public:
shared_ptr<PointList> vertices; //!< List of verticies
std::shared_ptr<PointList> vertices; //!< List of verticies
std::vector<int> A; //!< First vertex
std::vector<int> B; //!< Second vertex
std::vector<int> C; //!< Third vertex
@ -158,18 +158,18 @@ protected:
*/
struct MeshDataStruct {
std::string meshName;
shared_ptr<Mesh> mesh;
std::vector<shared_ptr<Variable> > vars;
std::shared_ptr<Mesh> mesh;
std::vector<std::shared_ptr<Variable> > vars;
};
//! Convert the mesh to a TriMesh (will return NULL if this is invalid)
shared_ptr<PointList> getPointList( shared_ptr<Mesh> mesh );
shared_ptr<TriMesh> getTriMesh( shared_ptr<Mesh> mesh );
shared_ptr<TriList> getTriList( shared_ptr<Mesh> mesh );
shared_ptr<const PointList> getPointList( shared_ptr<const Mesh> mesh );
shared_ptr<const TriMesh> getTriMesh( shared_ptr<const Mesh> mesh );
shared_ptr<const TriList> getTriList( shared_ptr<const Mesh> mesh );
std::shared_ptr<PointList> getPointList( std::shared_ptr<Mesh> mesh );
std::shared_ptr<TriMesh> getTriMesh( std::shared_ptr<Mesh> mesh );
std::shared_ptr<TriList> getTriList( std::shared_ptr<Mesh> mesh );
std::shared_ptr<const PointList> getPointList( std::shared_ptr<const Mesh> mesh );
std::shared_ptr<const TriMesh> getTriMesh( std::shared_ptr<const Mesh> mesh );
std::shared_ptr<const TriList> getTriList( std::shared_ptr<const Mesh> mesh );
} // IO namespace

7
IO/MeshDatabase.cpp
View File

@ -399,12 +399,13 @@ std::vector<MeshDatabase> read( const std::string& filename )
// Return the mesh type
IO::MeshType meshType( shared_ptr<IO::Mesh> mesh )
IO::MeshType meshType( std::shared_ptr<IO::Mesh> mesh )
{
IO::MeshType type = IO::Unknown;
if ( dynamic_pointer_cast<IO::PointList>(mesh).get()!=NULL ) {
if ( std::dynamic_pointer_cast<IO::PointList>(mesh).get()!=NULL ) {
type = IO::PointMesh;
} else if ( dynamic_pointer_cast<IO::TriList>(mesh).get()!=NULL || dynamic_pointer_cast<IO::TriMesh>(mesh).get()!=NULL ) {
} else if ( std::dynamic_pointer_cast<IO::TriList>(mesh).get()!=NULL ||
std::dynamic_pointer_cast<IO::TriMesh>(mesh).get()!=NULL ) {
type = IO::SurfaceMesh;
} 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( shared_ptr<IO::Mesh> mesh );
IO::MeshType meshType( std::shared_ptr<IO::Mesh> mesh );
} // IO namespace

14
IO/Reader.cpp
View File

@ -51,11 +51,11 @@ std::vector<IO::MeshDatabase> IO::getMeshList( const std::string& path, const st
// Read the given mesh domain
shared_ptr<IO::Mesh> IO::getMesh( const std::string& path, const std::string& timestep,
std::shared_ptr<IO::Mesh> IO::getMesh( const std::string& path, const std::string& timestep,
const IO::MeshDatabase& meshDatabase, int domain )
{
PROFILE_START("getMesh");
shared_ptr<IO::Mesh> mesh;
std::shared_ptr<IO::Mesh> mesh;
if ( meshDatabase.format==1 ) {
// Old format (binary doubles)
std::string filename = path + "/" + timestep + "/" + meshDatabase.domains[domain].file;
@ -72,7 +72,7 @@ shared_ptr<IO::Mesh> IO::getMesh( const std::string& path, const std::string& ti
ERROR("Error reading file");
if ( meshDatabase.type==IO::PointMesh ) {
size_t N = count/3;
shared_ptr<PointList> pointlist( new PointList(N) );
std::shared_ptr<PointList> pointlist( new PointList(N) );
std::vector<Point>& P = pointlist->points;
for (size_t i=0; i<N; i++) {
P[i].x = data[3*i+0];
@ -84,7 +84,7 @@ shared_ptr<IO::Mesh> IO::getMesh( const std::string& path, const std::string& ti
if ( count%9 != 0 )
ERROR("Error reading file (2)");
size_t N_tri = count/9;
shared_ptr<TriList> trilist( new TriList(N_tri) );
std::shared_ptr<TriList> trilist( new TriList(N_tri) );
std::vector<Point>& A = trilist->A;
std::vector<Point>& B = trilist->B;
std::vector<Point>& C = trilist->C;
@ -138,14 +138,14 @@ shared_ptr<IO::Mesh> IO::getMesh( const std::string& path, const std::string& ti
// Read the given variable for the given mesh domain
shared_ptr<IO::Variable> IO::getVariable( const std::string& path, const std::string& timestep,
std::shared_ptr<IO::Variable> IO::getVariable( const std::string& path, const std::string& timestep,
const MeshDatabase& meshDatabase, int domain, const std::string& variable )
{
std::pair<std::string,std::string> key(meshDatabase.domains[domain].name,variable);
std::map<std::pair<std::string,std::string>,DatabaseEntry>::const_iterator it;
it = meshDatabase.variable_data.find(key);
if ( it==meshDatabase.variable_data.end() )
return shared_ptr<IO::Variable>();
return std::shared_ptr<IO::Variable>();
const DatabaseEntry& database = it->second;
std::string filename = path + "/" + timestep + "/" + database.file;
FILE *fid = fopen(filename.c_str(),"rb");
@ -165,7 +165,7 @@ shared_ptr<IO::Variable> IO::getVariable( const std::string& path, const std::st
size_t count = fread(data,1,bytes,fid);
fclose(fid);
ASSERT(count==bytes);
shared_ptr<IO::Variable> var( new IO::Variable() );
std::shared_ptr<IO::Variable> var( new IO::Variable() );
var->dim = dim;
var->type = static_cast<IO::VariableType>(type);
var->name = variable;

4
IO/Reader.h
View File

@ -22,12 +22,12 @@ std::vector<IO::MeshDatabase> getMeshList( const std::string& path, const std::s
//! Read the given mesh domain
shared_ptr<IO::Mesh> getMesh( const std::string& path, const std::string& timestep,
std::shared_ptr<IO::Mesh> getMesh( const std::string& path, const std::string& timestep,
const MeshDatabase& meshDatabase, int domain );
//! Read the given mesh domain
shared_ptr<IO::Variable> getVariable( const std::string& path, const std::string& timestep,
std::shared_ptr<IO::Variable> getVariable( const std::string& path, const std::string& timestep,
const MeshDatabase& meshDatabase, int domain, const std::string& variable );

13
IO/Writer.cpp
View File

@ -26,7 +26,7 @@ static std::vector<IO::MeshDatabase> writeMeshesOrigFormat( const std::vector<IO
sprintf(fullpath,"%s/%s",path,filename);
FILE *fid = fopen(fullpath,"wb");
INSIST(fid!=NULL,std::string("Error opening file: ")+fullpath);
shared_ptr<IO::Mesh> mesh = meshData[i].mesh;
std::shared_ptr<IO::Mesh> mesh = meshData[i].mesh;
IO::MeshDatabase mesh_entry;
mesh_entry.name = meshData[i].meshName;
mesh_entry.type = meshType(mesh);
@ -42,18 +42,19 @@ 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 ( dynamic_pointer_cast<IO::PointList>(mesh).get()!=NULL ) {
if ( std::dynamic_pointer_cast<IO::PointList>(mesh).get()!=NULL ) {
// List of points
shared_ptr<IO::PointList> pointlist = dynamic_pointer_cast<IO::PointList>(mesh);
std::shared_ptr<IO::PointList> pointlist = std::dynamic_pointer_cast<IO::PointList>(mesh);
const std::vector<Point>& P = pointlist->points;
for (size_t i=0; i<P.size(); i++) {
double x[3];
x[0] = P[i].x; x[1] = P[i].y; x[2] = P[i].z;
fwrite(x,sizeof(double),3,fid);
}
} else if ( dynamic_pointer_cast<IO::TriList>(mesh).get()!=NULL || dynamic_pointer_cast<IO::TriMesh>(mesh).get()!=NULL ) {
} else if ( std::dynamic_pointer_cast<IO::TriList>(mesh).get()!=NULL ||
std::dynamic_pointer_cast<IO::TriMesh>(mesh).get()!=NULL ) {
// Triangle mesh
shared_ptr<IO::TriList> trilist = IO::getTriList(mesh);
std::shared_ptr<IO::TriList> trilist = IO::getTriList(mesh);
const std::vector<Point>& A = trilist->A;
const std::vector<Point>& B = trilist->B;
const std::vector<Point>& C = trilist->C;
@ -144,7 +145,7 @@ static std::vector<IO::MeshDatabase> writeMeshesNewFormat(
sprintf(fullpath,"%s/%s",path,filename);
FILE *fid = fopen(fullpath,"wb");
for (size_t i=0; i<meshData.size(); i++) {
shared_ptr<IO::Mesh> mesh = meshData[i].mesh;
std::shared_ptr<IO::Mesh> mesh = meshData[i].mesh;
meshes_written.push_back( write_domain(fid,filename,meshData[i],format) );
}
fclose(fid);

12
cmake/SharedPtr.cmake
View File

@ -58,6 +58,10 @@ FUNCTION( CONFIGURE_SHARED_PTR INSTALL_DIR NAMESPACE )
IF ( NOT NAMESPACE )
SET( NAMESPACE " " )
ENDIF()
SET( BOOST_SHARED_PTR 0 )
SET( MEMORY_SHARED_PTR 0 )
SET( MEMORY_TR1_SHARED_PTR 0 )
SET( TR1_MEMORY_TR1_SHARED_PTR 0 )
IF ( BOOST_SHARED_PTR )
FILE(WRITE "${CMAKE_CURRENT_BINARY_DIR}/tmp/shared_ptr.h" "#include \"boost/shared_ptr.hpp\"\n")
FILE(APPEND "${CMAKE_CURRENT_BINARY_DIR}/tmp/shared_ptr.h" "#include \"boost/weak_ptr.hpp\"\n")
@ -115,15 +119,18 @@ FUNCTION( WRITE_DUMMY_SHARED_PTR NAMESPACE FILENAME )
FILE(WRITE "${FILENAME}" "#ifndef DUMMY_SHARED_PTR_INC\n")
FILE(APPEND "${FILENAME}" "#define DUMMY_SHARED_PTR_INC\n")
FILE(APPEND "${FILENAME}" "namespace dummy {\n\n")
FILE(APPEND "${FILENAME}" "template<class T> void DefaultDeleter(T* p) {delete p;}\n\n")
FILE(APPEND "${FILENAME}" "template<class T> class shared_ptr {\n")
FILE(APPEND "${FILENAME}" "public:\n")
FILE(APPEND "${FILENAME}" " shared_ptr( ): obj(NULL), count(NULL) {}\n")
FILE(APPEND "${FILENAME}" " shared_ptr( T *ptr ): obj(ptr), count(NULL) { if (ptr) { count = new int; (*count)=1; } } \n")
FILE(APPEND "${FILENAME}" " shared_ptr( T *ptr, void (*D)(T*)=DefaultDeleter<T>):\n")
FILE(APPEND "${FILENAME}" " obj(ptr), deleter(D), count(NULL) { if (ptr) { count = new int; (*count)=1; } } \n")
FILE(APPEND "${FILENAME}" " shared_ptr( const shared_ptr<T>& rhs ): \n")
FILE(APPEND "${FILENAME}" " obj(rhs.get()), count(rhs.count) { if ( count!=NULL ) { ++(*count); } } \n")
FILE(APPEND "${FILENAME}" " template<class U> shared_ptr( const shared_ptr<U>& rhs ): \n")
FILE(APPEND "${FILENAME}" " obj(rhs.get()), count(rhs.count) { if ( count!=NULL ) { ++(*count); } } \n")
FILE(APPEND "${FILENAME}" " shared_ptr& operator=( const shared_ptr<T>& rhs ) { obj=rhs.obj; count=rhs.count; ++(*count); return *this; } \n")
FILE(APPEND "${FILENAME}" " shared_ptr& operator=( const shared_ptr<T>& rhs )\n")
FILE(APPEND "${FILENAME}" " { obj=rhs.obj; count=rhs.count; ++(*count); return *this; } \n")
FILE(APPEND "${FILENAME}" " ~shared_ptr( ) { reset(); }\n")
FILE(APPEND "${FILENAME}" " void reset( T *ptr ) { reset(); obj=ptr; count=new int; (*count)=1; }\n")
FILE(APPEND "${FILENAME}" " void reset( void ) { \n")
@ -137,6 +144,7 @@ FUNCTION( WRITE_DUMMY_SHARED_PTR NAMESPACE FILENAME )
FILE(APPEND "${FILENAME}" " bool operator!=( const T * rhs ) const { return obj!=rhs; } \n")
FILE(APPEND "${FILENAME}" "protected:\n")
FILE(APPEND "${FILENAME}" " T *obj;\n")
FILE(APPEND "${FILENAME}" " void (*deleter)(T*);\n")
FILE(APPEND "${FILENAME}" " volatile int *count;\n")
FILE(APPEND "${FILENAME}" "template<class T1, class U> friend shared_ptr<T1> dynamic_pointer_cast( shared_ptr<U> const & );\n")
FILE(APPEND "${FILENAME}" "template<class T1, class U> friend shared_ptr<T1> const_pointer_cast( shared_ptr<U> const & );\n")

258
common/Array.cpp
View File

@ -1,258 +0,0 @@
#if 0
#include "Array.h"
/*
* Array.cpp
*
* Created by James Mcclure on 3/31/09.
* Copyright 2009 __MyCompanyName__. All rights reserved.
*
*/
// *****************************************
// ******** class IntArray ****************
// *****************************************
IntArray::IntArray()
{
m=n=o=Length=0;
}
IntArray::IntArray(int size)
{
data = new int [size];
Length = size;
m = size;
n = 1;
o = 1;
}
IntArray::IntArray(int nx, int ny)
{
m = nx;
n = ny;
o = 1;
Length = m*n;
data = new int [Length];
}
IntArray::IntArray(int nx, int ny, int nz)
{
m = nx;
n = ny;
o = nz;
Length = m*n*o;
data = new int [Length];
}
IntArray::~IntArray()
{
delete data;
}
void IntArray::New(int size)
{
m=size;
n = 1;
o = 1;
data = new int [size];
Length = size;
}
void IntArray::New(int nx, int ny)
{
m = nx;
n = ny;
o = 1;
Length = m*n;
data = new int [Length];
}
void IntArray::New(int nx, int ny,int nz)
{
m = nx;
n = ny;
o = nz;
Length = m*n*o;
data = new int [Length];
}
int IntArray::e(int i)
{
return data[i];
}
int IntArray::e(int i, int j)
{
return data[m*j+i];
}
int IntArray::e(int i, int j, int k)
{
return data[m*n*k+m*j+i];
}
// *****************************************
// ******** class DoubleArray **************
// *****************************************
DoubleArray::DoubleArray()
{
m=n=o=Length=0;
}
DoubleArray::DoubleArray(int size)
{
m=size;
n = 1;
o = 1;
data = new double [size];
Length = size;
}
DoubleArray::DoubleArray(int nx, int ny)
{
m = nx;
n = ny;
o = 1;
Length = m*n;
data = new double [Length];
}
DoubleArray::DoubleArray(int nx, int ny,int nz)
{
m = nx;
n = ny;
o = nz;
Length = m*n*o;
data = new double [Length];
}
void DoubleArray::New(int size)
{
m=size;
n = 1;
o = 1;
data = new double [size];
Length = size;
}
void DoubleArray::New(int nx, int ny)
{
m = nx;
n = ny;
o = 1;
Length = m*n;
data = new double [Length];
}
void DoubleArray::New(int nx, int ny,int nz)
{
m = nx;
n = ny;
o = nz;
Length = m*n*o;
data = new double [Length];
}
DoubleArray::~DoubleArray()
{
delete data;
}
double DoubleArray::e(int i)
{
return data[i];
}
double DoubleArray::e(int i, int j)
{
return data[j*m+i];
}
double DoubleArray::e(int i, int j, int k)
{
return data[k*m*n+j*m+i];
}
extern DoubleArray IncreaseSize(DoubleArray &A, int addLength)
{
if (addLength<0) {
printf("IncreaseSize(Array,Length)","Length needs to be >0.");
return DoubleArray();
}
int newM,newN,newO;
if (A.o>1) {
if (addLength%(A.m*A.n)!=0) {
printf("IncreaseSize(Array,Length)","Length needs to be a multiple of m*n");
return DoubleArray();
}
newM = A.m;
newN = A.n;
newO = A.o + addLength/(A.m*A.n);
}
else if (A.n>1) {
if (addLength%(A.m)!=0) {
printf("IncreaseSize(Array,Length)","Length needs to be a multiple of m");
return DoubleArray();
}
newM = A.m;
newN = A.n + addLength/A.m;
newO = 1;
}
else {
newM = A.m + addLength;
newN = 1;
newO = 1;
}
DoubleArray toReturn(newM,newN,newO);
memcpy(toReturn.Pointer(),A.Pointer(),A.Length*sizeof(double));
return toReturn;
}
extern IntArray IncreaseSize(IntArray &A, int addLength)
{
if (addLength<0) {
printf("IncreaseSize(Array,Length)","Length needs to be >0.");
return IntArray();
}
int newM,newN,newO;
if (A.o>1) {
if (addLength%(A.m*A.n)!=0) {
printf("IncreaseSize(Array,Length)","Length needs to be a multiple of m*n");
return IntArray();
}
newM = A.m;
newN = A.n;
newO = A.o + addLength/(A.m*A.n);
}
else if (A.n>1) {
if (addLength%(A.m)!=0) {
printf("IncreaseSize(Array,Length)","Length needs to be a multiple of m");
return IntArray();
}
newM = A.m;
newN = A.n + addLength/A.m;
newO = 1;
}
else {
newM = A.m + addLength;
newN = 1;
newO = 1;
}
IntArray toReturn(newM,newN,newO);
memcpy(toReturn.Pointer(),A.Pointer(),A.Length*sizeof(int));
return toReturn;
}
DoubleArray DoubleArray::Copy()
{
DoubleArray CopyInto(m,n,o);
// Check that the allocation worked.
if (CopyInto.Length!=Length) return CopyInto; // Failed. Already printed an error message.
memcpy(CopyInto.Pointer(),Pointer(),Length*sizeof(double));
return CopyInto;
}
#endif

705
common/Array.h
View File

@ -1,413 +1,330 @@
#ifndef ARRAY_H_INC
#define ARRAY_H_INC
#ifndef included_ArrayClass
#define included_ArrayClass
#include <iostream>
#include <string.h>
#include <memory>
#include <vector>
#include "shared_ptr.h"
#include "common/Utilities.h"
// Define a macro to check if the index is valid
// Only perform the check if we are compiling in debug mode
#define GET_ARRAY_INDEX(i1,i2,i3) i1+d_N[0]*(i2+d_N[1]*i3)
#ifdef DEBUG
#ifdef USE_CUDA
#define CHECK_INDEX(i,j,k) \
if ( (i+j*m+k*m*n)<0 || (i+j*m+k*m*n)>=Length ) { \
printf("Index is out of bounds\n"); }
#else
#include "common/Utilities.h"
#define CHECK_INDEX(i,j,k) \
if ( (i+j*m+k*m*n)<0 || (i+j*m+k*m*n)>=Length ) { \
ERROR("Index is out of bounds\n"); }
#endif
#define CHECK_ARRAY_INDEX(i1,i2,i3) \
if ( GET_ARRAY_INDEX(i1,i2,i3)<0 || GET_ARRAY_INDEX(i1,i2,i3)>d_length ) \
ERROR("Index exceeds array bounds");
#else
#define CHECK_INDEX(i,j,k)
#define CHECK_ARRAY_INDEX(i1,i2,i3)
#endif
// ********** ARRAY CLASS INFO **************************************
/*
//..............................................................
// Overview of array classes in Array.h
//......... Declaration of array objects........................
// supports integer data (IntArray) or double data (DoubleArray)
// supports one, two or three dimensional arrays
IntArray A(m); // array of integers, Length n
IntArray A(m,n); // array of integers, dimensions: m x n
DoubleArray A(m,n,o); // array of doubles, dimensions m x n x o
//............ Access the size of the array.....................
A.m; // size of first dimension
A.n; // size of second dimension
A.o; // size of third dimension
A.Length; // total number of values stored
//........... Access the array entries .........................
A(i); // return data[i]
A(i,j); // return data[j*m+i]
A(i,j,k); // return data[k*m*n+j*m+i]
//..............................................................
/*!
* Class Array is a simple array class
*/
using namespace std;
class IntArray{
template<class TYPE>
class Array
{
public:
IntArray Copy();
int Length;
int m,n,o;
int *data;
IntArray();
IntArray(int size);
IntArray(int nx,int ny);
IntArray(int nx,int ny,int nz);
~IntArray();
void New(int size);
void New(int nx, int ny);
void New(int nx, int ny, int nz);
int & operator()(int i)
{ CHECK_INDEX(i,0,0) return data[i];}
int & operator()(int i, int j)
{ CHECK_INDEX(i,j,0) return data[j*m+i];}
int & operator()(int i, int j, int k)
{ CHECK_INDEX(i,j,k) return data[k*m*n+j*m+i];}
int e(int i);
int e(int i, int j);
int e(int i, int j, int k);
int *Pointer() {return data;}
/*!
* Create a new empty Array
*/
Array( );
/*!
* Create a new 1D Array with the given number of elements
* @param N Number of elements in the array
*/
Array( size_t N );
/*!
* Create a new 2D Array with the given number of rows and columns
* @param N_rows Number of rows
* @param N_columns Number of columns
*/
Array( size_t N_rows, size_t N_columns );
/*!
* Create a new 3D Array with the given number of rows and columns
* @param N1 Number of rows
* @param N2 Number of columns
* @param N3 Number of elements in the third dimension
*/
Array( size_t N1, size_t N2, size_t N3 );
/*!
* Create a multi-dimensional Array with the given number of elements
* @param N Number of elements in each dimension
*/
Array( const std::vector<size_t>& N );
/*!
* Copy constructor
* @param rhs Array to copy
*/
Array( const Array& rhs );
/*!
* Assignment operator
* @param rhs Array to copy
*/
Array& operator=( const Array& rhs );
/*!
* Create a 1D Array view to a raw block of data
* @param N Number of elements in the array
* @param data Pointer to the data
*/
static std::shared_ptr<Array> view( size_t N, std::shared_ptr<TYPE> data );
/*!
* Create a new 2D Array with the given number of rows and columns
* @param N_rows Number of rows
* @param N_columns Number of columns
* @param data Pointer to the data
*/
static std::shared_ptr<Array> view( size_t N_rows, size_t N_columns, std::shared_ptr<TYPE> data );
/*!
* Create a new 3D Array view to a raw block of data
* @param N1 Number of rows
* @param N2 Number of columns
* @param N3 Number of elements in the third dimension
* @param data Pointer to the data
*/
static std::shared_ptr<Array> view( size_t N1, size_t N2, size_t N3, std::shared_ptr<TYPE> data );
/*!
* Create a multi-dimensional Array view to a raw block of data
* @param N Number of elements in each dimension
*/
static std::shared_ptr<Array> view( const std::vector<size_t>& N, std::shared_ptr<TYPE> data );
/*!
* Create a 1D Array view to a raw block of data
* @param N Number of elements in the array
* @param data Pointer to the data
*/
static std::shared_ptr<const Array> constView( size_t N, std::shared_ptr<const TYPE> data );
/*!
* Create a new 2D Array with the given number of rows and columns
* @param N_rows Number of rows
* @param N_columns Number of columns
* @param data Pointer to the data
*/
static std::shared_ptr<const Array> constView( size_t N_rows, size_t N_columns, std::shared_ptr<const TYPE> data );
/*!
* Create a new 3D Array view to a raw block of data
* @param N1 Number of rows
* @param N2 Number of columns
* @param N3 Number of elements in the third dimension
* @param data Pointer to the data
*/
static std::shared_ptr<const Array> constView( size_t N1, size_t N2, size_t N3, std::shared_ptr<const TYPE> data );
/*!
* Create a multi-dimensional Array view to a raw block of data
* @param N Number of elements in each dimension
*/
static std::shared_ptr<const Array> constView( const std::vector<size_t>& N, std::shared_ptr<const TYPE> data );
/*!
* Convert an array of one type to another. This may or may not allocate new memory.
* @param array Input array
*/
template<class TYPE2>
static std::shared_ptr<Array<TYPE2> > convert( std::shared_ptr<Array<TYPE> > array );
/*!
* Convert an array of one type to another. This may or may not allocate new memory.
* @param array Input array
*/
template<class TYPE2>
static std::shared_ptr<const Array<TYPE2> > convert( std::shared_ptr<const Array<TYPE> > array );
/*!
* Copy and convert data from another array to this array
* @param array Source array
*/
template<class TYPE2>
void copy( const Array<TYPE2>& array );
/*!
* Copy and convert data from a raw vector to this array.
* Note: The current array must be allocated to the proper size first.
* @param array Source array
*/
template<class TYPE2>
void copy( const TYPE2* array );
/*!
* Fill the array with the given value
* @param value Value to fill
*/
void fill( const TYPE& value );
//! Destructor
~Array( );
//! Return the size of the Array
inline int ndim( ) const { return d_ndim; }
//! Return the size of the Array
inline std::vector<size_t> size( ) const { return std::vector<size_t>(d_N,d_N+d_ndim); }
//! Return the size of the Array
inline size_t size( int d ) const { return d_N[d]; }
//! Return the size of the Array
inline size_t length( ) const { return d_length; }
//! Return true if the Array is empty
inline bool empty( ) const { return d_length==0; }
/*!
* Resize the Array
* @param N NUmber of elements
*/
void resize( size_t N );
/*!
* Resize the Array
* @param N_rows Number of rows
* @param N_columns Number of columns
*/
void resize( size_t N_rows, size_t N_columns );
/*!
* Resize the Array
* @param N1 Number of rows
* @param N2 Number of columns
* @param N3 Number of elements in the third dimension
*/
void resize( size_t N1, size_t N2, size_t N3 );
/*!
* Resize the Array
* @param N Number of elements in each dimension
*/
void resize( const std::vector<size_t>& N );
/*!
* Reshape the Array (total size of array will not change)
* @param N Number of elements in each dimension
*/
void reshape( const std::vector<size_t>& N );
/*!
* Access the desired element
* @param i The row index
* @param j The column index
*/
inline TYPE& operator()( size_t i ) { CHECK_ARRAY_INDEX(i,0,0) return d_data[i]; }
/*!
* Access the desired element
* @param i The row index
* @param j The column index
*/
inline const TYPE& operator()( size_t i ) const { CHECK_ARRAY_INDEX(i,0,0) return d_data[i]; }
/*!
* Access the desired element
* @param i The row index
* @param j The column index
*/
inline TYPE& operator()( size_t i, size_t j ) { CHECK_ARRAY_INDEX(i,j,0) return d_data[i+j*d_N[0]]; }
/*!
* Access the desired element
* @param i The row index
* @param j The column index
*/
inline const TYPE& operator()( size_t i, size_t j ) const { CHECK_ARRAY_INDEX(i,j,0) return d_data[i+j*d_N[0]]; }
/*!
* Access the desired element
* @param i The row index
* @param j The column index
*/
inline TYPE& operator()( size_t i, size_t j, size_t k ) { CHECK_ARRAY_INDEX(i,j,k) return d_data[GET_ARRAY_INDEX(i,j,k)]; }
/*!
* Access the desired element
* @param i The row index
* @param j The column index
*/
inline const TYPE& operator()( size_t i, size_t j, size_t k ) const { CHECK_ARRAY_INDEX(i,j,k) return d_data[GET_ARRAY_INDEX(i,j,k)]; }
//! Check if two matricies are equal
bool operator==( const Array& rhs ) const;
//! Check if two matricies are not equal
inline bool operator!=( const Array& rhs ) const { return !this->operator==(rhs); }
//! Return the pointer to the raw data
inline std::shared_ptr<TYPE> getPtr( ) { return d_ptr; }
//! Return the pointer to the raw data
inline std::shared_ptr<const TYPE> getPtr( ) const { return d_ptr; }
//! Return the pointer to the raw data
inline TYPE* get( ) { return d_data; }
//! Return the pointer to the raw data
inline const TYPE* get( ) const { return d_data; }
//! Return true if NaNs are present
inline bool NaNs( ) const;
//! Return the smallest value
inline TYPE min( ) const;
//! Return the largest value
inline TYPE max( ) const;
//! Return the sum of all elements
inline TYPE sum( ) const;
//! Return the sum of all elements in a given direction
std::shared_ptr<Array<TYPE> > sum( int dir ) const;
private:
int d_ndim;
size_t d_N[3];
size_t d_length;
TYPE *d_data;
std::shared_ptr<TYPE> d_ptr;
void allocate( const std::vector<size_t>& N );
};
class DoubleArray{
public:
DoubleArray Copy();
int Length;
int m;
int n;
int o;
double *data;
DoubleArray();
DoubleArray(int size);
DoubleArray(int nx, int ny);
DoubleArray(int nx, int ny, int nz);
~DoubleArray();
void New(int size);
void New(int nx, int ny);
void New(int nx, int ny, int nz);
double & operator()(int i)
{ CHECK_INDEX(i,0,0) return data[i];}
double & operator()(int i, int j)
{ CHECK_INDEX(i,j,0) return data[j*m+i];}
double & operator()(int i, int j, int k)
{ CHECK_INDEX(i,j,k) return data[k*m*n+j*m+i];}
double e(int i);
double e(int i, int j);
double e(int i, int j, int k);
double *Pointer() {return data;}
};
/*
* Array.cpp
*
* Created by James Mcclure on 3/31/09.
* Copyright 2009 __MyCompanyName__. All rights reserved.
*
*/
// *****************************************
// ******** class IntArray ****************
// *****************************************
/*IntArray::IntArray();
IntArray::IntArray(int size);
IntArray::IntArray(int nx, int ny);
IntArray::IntArray(int nx, int ny, int nz);
IntArray::~IntArray();
void IntArray::New(int size);
void IntArray::New(int nx, int ny);
void IntArray::New(int nx, int ny,int nz);
int IntArray::e(int i);
int IntArray::e(int i, int j);
int IntArray::e(int i, int j, int k);
// *****************************************
// ******** class DoubleArray **************
// *****************************************
DoubleArray::DoubleArray();
DoubleArray::DoubleArray(int size);
DoubleArray::DoubleArray(int nx, int ny);
DoubleArray::DoubleArray(int nx, int ny,int nz);
DoubleArray::~DoubleArray();
void DoubleArray::New(int size);
void DoubleArray::New(int nx, int ny);
void DoubleArray::New(int nx, int ny,int nz);
double DoubleArray::e(int i);
double DoubleArray::e(int i, int j);
double DoubleArray::e(int i, int j, int k);
extern DoubleArray IncreaseSize(DoubleArray &A, int addLength);
extern IntArray IncreaseSize(IntArray &A, int addLength);
*/
/*
* Array.cpp
*
* Created by James Mcclure on 3/31/09.
* Copyright 2009 __MyCompanyName__. All rights reserved.
*
*/
/*
* Array.cpp
*
* Created by James Mcclure on 3/31/09.
* Copyright 2009 __MyCompanyName__. All rights reserved.
*
*/
IntArray::IntArray()
{
m=n=o=Length=0;
}
IntArray::IntArray(int size)
{
data = new int [size];
Length = size;
m = size;
n = 1;
o = 1;
}
IntArray::IntArray(int nx, int ny)
{
m = nx;
n = ny;
o = 1;
Length = m*n;
data = new int [Length];
}
IntArray::IntArray(int nx, int ny, int nz)
{
m = nx;
n = ny;
o = nz;
Length = m*n*o;
data = new int [Length];
}
IntArray::~IntArray()
{
delete[] data;
}
typedef Array<int> IntArray;
typedef Array<double> DoubleArray;
void IntArray::New(int size)
{
m=size;
n = 1;
o = 1;
data = new int [size];
Length = size;
}
void IntArray::New(int nx, int ny)
{
m = nx;
n = ny;
o = 1;
Length = m*n;
data = new int [Length];
}
void IntArray::New(int nx, int ny,int nz)
{
m = nx;
n = ny;
o = nz;
Length = m*n*o;
data = new int [Length];
}
int IntArray::e(int i)
{
return data[i];
}
int IntArray::e(int i, int j)
{
return data[m*j+i];
}
int IntArray::e(int i, int j, int k)
{
return data[m*n*k+m*j+i];
}
// *****************************************
// ******** class DoubleArray **************
// *****************************************
DoubleArray::DoubleArray()
{
m=n=o=Length=0;
}
DoubleArray::DoubleArray(int size)
{
m=size;
n = 1;
o = 1;
data = new double [size];
Length = size;
}
DoubleArray::DoubleArray(int nx, int ny)
{
m = nx;
n = ny;
o = 1;
Length = m*n;
data = new double [Length];
}
DoubleArray::DoubleArray(int nx, int ny,int nz)
{
m = nx;
n = ny;
o = nz;
Length = m*n*o;
data = new double [Length];
}
void DoubleArray::New(int size)
{
m=size;
n = 1;
o = 1;
data = new double [size];
Length = size;
}
void DoubleArray::New(int nx, int ny)
{
m = nx;
n = ny;
o = 1;
Length = m*n;
data = new double [Length];
}
void DoubleArray::New(int nx, int ny,int nz)
{
m = nx;
n = ny;
o = nz;
Length = m*n*o;
data = new double [Length];
}
DoubleArray::~DoubleArray()
{
delete [] data;
}
double DoubleArray::e(int i)
{
return data[i];
}
double DoubleArray::e(int i, int j)
{
return data[j*m+i];
}
double DoubleArray::e(int i, int j, int k)
{
return data[k*m*n+j*m+i];
}
extern DoubleArray IncreaseSize(DoubleArray &A, int addLength)
{
if (addLength<0) {
printf("IncreaseSize(Array,Length) Length needs to be >0.");
return DoubleArray();
}
int newM,newN,newO;
if (A.o>1) {
if (addLength%(A.m*A.n)!=0) {
printf("IncreaseSize(Array,Length) Length needs to be a multiple of m*n");
return DoubleArray();
}
newM = A.m;
newN = A.n;
newO = A.o + addLength/(A.m*A.n);
}
else if (A.n>1) {
if (addLength%(A.m)!=0) {
printf("IncreaseSize(Array,Length) Length needs to be a multiple of m");
return DoubleArray();
}
newM = A.m;
newN = A.n + addLength/A.m;
newO = 1;
}
else {
newM = A.m + addLength;
newN = 1;
newO = 1;
}
DoubleArray toReturn(newM,newN,newO);
memcpy(toReturn.Pointer(),A.Pointer(),A.Length*sizeof(double));
return toReturn;
}
extern IntArray IncreaseSize(IntArray &A, int addLength)
{
if (addLength<0) {
printf("IncreaseSize(Array,Length) Length needs to be >0.");
return IntArray();
}
int newM,newN,newO;
if (A.o>1) {
if (addLength%(A.m*A.n)!=0) {
printf("IncreaseSize(Array,Length) Length needs to be a multiple of m*n");
return IntArray();
}
newM = A.m;
newN = A.n;
newO = A.o + addLength/(A.m*A.n);
}
else if (A.n>1) {
if (addLength%(A.m)!=0) {
printf("IncreaseSize(Array,Length) Length needs to be a multiple of m");
return IntArray();
}
newM = A.m;
newN = A.n + addLength/A.m;
newO = 1;
}
else {
newM = A.m + addLength;
newN = 1;
newO = 1;
}
IntArray toReturn(newM,newN,newO);
memcpy(toReturn.Pointer(),A.Pointer(),A.Length*sizeof(int));
return toReturn;
}
DoubleArray DoubleArray::Copy()
{
DoubleArray CopyInto(m,n,o);
// Check that the allocation worked.
if (CopyInto.Length!=Length) return CopyInto; // Failed. Already printed an error message.
memcpy(CopyInto.Pointer(),Pointer(),Length*sizeof(double));
return CopyInto;
}
#include "common/Array.hpp"
#endif

377
common/Array.hpp Normal file
View File

@ -0,0 +1,377 @@
#ifndef included_ArrayClass_hpp
#define included_ArrayClass_hpp
#include "common/Array.h"
#include "common/Utilities.h"
#include <algorithm>
#include <limits>
template<class TYPE>
void DeleteArray( TYPE* x )
{
delete [] x;
}
/********************************************************
* Constructors *
********************************************************/
template<class TYPE>
Array<TYPE>::Array( )
{
d_ndim = 0;
d_length = 0;
for (size_t i=0; i<sizeof(d_N)/sizeof(size_t); i++)
d_N[i] = 1;
d_N[0] = 0;
d_data = d_ptr.get();
}
template<class TYPE>
Array<TYPE>::Array( size_t N )
{
allocate(std::vector<size_t>(1,N));
}
template<class TYPE>
Array<TYPE>::Array( size_t N_rows, size_t N_columns )
{
std::vector<size_t> N(2);
N[0] = N_rows;
N[1] = N_columns;
allocate(N);
}
template<class TYPE>
Array<TYPE>::Array( size_t N1, size_t N2, size_t N3 )
{
std::vector<size_t> N(3);
N[0] = N1;
N[1] = N2;
N[2] = N3;
allocate(N);
}
template<class TYPE>
Array<TYPE>::Array( const std::vector<size_t>& N )
{
allocate(N);
}
template<class TYPE>
void Array<TYPE>::allocate( const std::vector<size_t>& N )
{
d_ndim = N.size();
d_length = 1;
for (size_t i=0; i<sizeof(d_N)/sizeof(size_t); i++)
d_N[i] = 1;
for (size_t i=0; i<N.size(); i++) {
d_N[i] = N[i];
d_length *= N[i];
}
if ( N.empty() ) {
d_N[0] = 0;
d_length = 0;
}
if ( d_length==0 )
d_ptr.reset();
else
d_ptr = std::shared_ptr<TYPE>(new TYPE[d_length],DeleteArray<TYPE>);
d_data = d_ptr.get();
if ( d_length>0 && d_data==NULL )
ERROR("Failed to allocate array");
}
template<class TYPE>
Array<TYPE>::Array( const Array& rhs ):
d_ndim(rhs.d_ndim), d_length(rhs.d_length), d_data(NULL)
{
allocate( std::vector<size_t>(rhs.d_N,rhs.d_N+rhs.d_ndim) );
for (size_t i=0; i<d_length; i++)
d_data[i] = rhs.d_data[i];
}
template<class TYPE>
Array<TYPE>& Array<TYPE>::operator=( const Array& rhs )
{
if ( this == &rhs )
return *this;
this->allocate( std::vector<size_t>(rhs.d_N,rhs.d_N+rhs.d_data) );
for (size_t i=0; i<d_length; i++)
this->d_data[i] = rhs.d_data[i];
return *this;
}
template<class TYPE>
Array<TYPE>::~Array( )
{
}
/********************************************************
* Resize the array *
********************************************************/
template<class TYPE>
void Array<TYPE>::resize( size_t N )
{
resize(std::vector<size_t>(1,N));
}
template<class TYPE>
void Array<TYPE>::resize( size_t N1, size_t N2 )
{
std::vector<size_t> N(2);
N[0] = N1;
N[1] = N2;
resize(N);
}
template<class TYPE>
void Array<TYPE>::resize( size_t N1, size_t N2, size_t N3 )
{
std::vector<size_t> N(3);
N[0] = N1;
N[1] = N2;
N[2] = N3;
resize(N);
}
template<class TYPE>
void Array<TYPE>::resize( const std::vector<size_t>& N )
{
// Check if the array actually changed size
size_t new_length = 1;
for (size_t i=0; i<N.size(); i++)
new_length *= N[i];
bool changed = new_length!=d_length;
for (size_t i=0; i<N.size(); i++)
changed = changed || N[i]!=d_N[i];
if ( !changed )
return;
// Store the old data
const size_t ndim_max = sizeof(d_N)/sizeof(size_t);
std::vector<size_t> N1(ndim_max,1), N2(ndim_max,1);
for (size_t d=0; d<d_ndim; d++)
N1[d] = d_N[d];
for (size_t d=0; d<N.size(); d++)
N2[d] = N[d];
if ( d_ndim==0 ) { N1[0] = 0; }
if ( N.empty() ) { N2[0] = 0; }
std::shared_ptr<TYPE> old_data = d_ptr;
// Allocate new data
allocate(N);
// Copy the old values
if ( d_length > 0 ) {
ASSERT(sizeof(d_N)/sizeof(size_t)==3);
TYPE *data1 = old_data.get();
TYPE *data2 = d_data;
for (size_t k=0; k<std::min(N1[2],N2[2]); k++) {
for (size_t j=0; j<std::min(N1[1],N2[1]); j++) {
for (size_t i=0; i<std::min(N1[0],N2[0]); i++) {
size_t index1 = i + j*N1[0] + k*N1[0]*N1[1];
size_t index2 = i + j*N2[0] + k*N2[0]*N2[1];
data2[index2] = data1[index1];
}
}
}
}
}
/********************************************************
* Rehape the array *
********************************************************/
template<class TYPE>
void Array<TYPE>::reshape( const std::vector<size_t>& N )
{
size_t new_length = 1;
for (size_t i=0; i<N.size(); i++)
new_length *= N[i];
if ( new_length!=d_length )
ERROR("reshape is not allowed to change the array size");
d_ndim = N.size();
for (size_t i=0; i<sizeof(d_N)/sizeof(size_t); i++)
d_N[i] = 1;
for (size_t i=0; i<N.size(); i++)
d_N[i] = N[i];
}
/********************************************************
* Operator overloading *
********************************************************/
template<class TYPE>
bool Array<TYPE>::operator==( const Array& rhs ) const
{
if ( this==&rhs )
return true;
if ( d_length!=rhs.d_length )
return false;
bool match = true;
for (size_t i=0; i<d_length; i++)
match = match && d_data[i]==rhs.d_data[i];
return match;
}
/********************************************************
* Get a view of an C array *
********************************************************/
template<class TYPE>
std::shared_ptr<Array<TYPE> > Array<TYPE>::view( size_t N, std::shared_ptr<TYPE> data )
{
view(std::vector<size_t>(1,N),data);
}
template<class TYPE>
std::shared_ptr<Array<TYPE> > Array<TYPE>::view( size_t N1, size_t N2, std::shared_ptr<TYPE> data )
{
std::vector<size_t> N(2);
N[0] = N1;
N[1] = N2;
view(N,data);
}
template<class TYPE>
std::shared_ptr<Array<TYPE> > Array<TYPE>::view( size_t N1, size_t N2, size_t N3, std::shared_ptr<TYPE> data )
{
std::vector<size_t> N(3);
N[0] = N1;
N[1] = N2;
N[2] = N3;
view(N,data);
}
template<class TYPE>
std::shared_ptr<const Array<TYPE> > Array<TYPE>::constView( size_t N, std::shared_ptr<const TYPE> data )
{
constView(std::vector<size_t>(1,N),data);
}
template<class TYPE>
std::shared_ptr<const Array<TYPE> > Array<TYPE>::constView( size_t N1, size_t N2, std::shared_ptr<const TYPE> data )
{
std::vector<size_t> N(2);
N[0] = N1;
N[1] = N2;
constView(N,data);
}
template<class TYPE>
std::shared_ptr<const Array<TYPE> > Array<TYPE>::constView( size_t N1, size_t N2, size_t N3, std::shared_ptr<const TYPE> data )
{
std::vector<size_t> N(3);
N[0] = N1;
N[1] = N2;
N[2] = N3;
constView(N,data);
}
template<class TYPE>
std::shared_ptr<Array<TYPE> > Array<TYPE>::view( const std::vector<size_t>& N, std::shared_ptr<TYPE> data )
{
std::shared_ptr<Array<TYPE> > array(new Array<TYPE>());
array->d_ndim = N.size();
array->d_length = 1;
for (size_t i=0; i<N.size(); i++) {
array->d_N[i] = N[i];
array->d_length *= N[i];
}
array->d_ptr = data;
array->d_data = array->d_ptr.get();
return array;
}
template<class TYPE>
std::shared_ptr<const Array<TYPE> > Array<TYPE>::constView( const std::vector<size_t>& N, std::shared_ptr<const TYPE> data )
{
return view(N,std::const_pointer_cast<TYPE>(data));
}
/********************************************************
* Convert array types *
********************************************************/
template<class TYPE>
template<class TYPE2>
std::shared_ptr<Array<TYPE2> > Array<TYPE>::convert( std::shared_ptr<Array<TYPE> > array )
{
std::shared_ptr<Array<TYPE2> > array2( new Array<TYPE2>(array->size()) );
array2.copy( *array );
return array2;
}
template<class TYPE>
template<class TYPE2>
std::shared_ptr<const Array<TYPE2> > Array<TYPE>::convert( std::shared_ptr<const Array<TYPE> > array )
{
return Array<TYPE>::convert( std::const_pointer_cast<Array<TYPE2> >(array) );
}
template<class TYPE>
template<class TYPE2>
void Array<TYPE>::copy( const Array<TYPE2>& array )
{
resize( std::vector<size_t>(array.d_N,array.d_N+array.d_ndim) );
const TYPE2 *src = array.d_data;
for (size_t i=0; i<d_length; i++)
d_data[i] = static_cast<TYPE>(src[i]);
}
template<class TYPE>
template<class TYPE2>
void Array<TYPE>::copy( const TYPE2* src )
{
for (size_t i=0; i<d_length; i++)
d_data[i] = static_cast<TYPE>(src[i]);
}
template<class TYPE>
void Array<TYPE>::fill( const TYPE& value )
{
for (size_t i=0; i<d_length; i++)
d_data[i] = value;
}
/********************************************************
* Simple math operations *
********************************************************/
template<class TYPE>
bool Array<TYPE>::NaNs( ) const
{
bool test = false;
for (size_t i=0; i<d_length; i++)
test = test || d_data[i]!=d_data[i];
return test;
}
template<class TYPE>
TYPE Array<TYPE>::min( ) const
{
TYPE x = std::numeric_limits<TYPE>::max();
for (size_t i=0; i<d_length; i++)
x = std::min(x,d_data[i]);
return x;
}
template<class TYPE>
TYPE Array<TYPE>::max( ) const
{
TYPE x = std::numeric_limits<TYPE>::min();
for (size_t i=0; i<d_length; i++)
x = std::max(x,d_data[i]);
return x;
}
template<class TYPE>
TYPE Array<TYPE>::sum( ) const
{
TYPE x = 0;
for (size_t i=0; i<d_length; i++)
x += d_data[i];
return x;
}
template<class TYPE>
std::shared_ptr<Array<TYPE> > Array<TYPE>::sum( int dir ) const
{
std::vector<size_t> size_ans = size();
std::shared_ptr<Array<TYPE> > ans( new Array<TYPE>(size_ans) );
size_t N1=1, N2=1, N3=1;
for (int d=0; d<std::min(dir,d_ndim); d++)
N1 *= d_N[d];
N2 = d_N[dir];
for (int d=dir+1; d<std::min(dir,d_ndim); d++)
N3 *= d_N[d];
TYPE* data2 = ans->d_data;
for (int i3=0; i3<N3; i3++) {
for (int i1=0; i1<N1; i1++) {
TYPE x = 0;
for (size_t i2=0; i2<N2; i2++)
x += d_data[i1+i2*N1+i3*N1*N2];
data2[i1+i3*N1] = x;
}
}
return ans;
}
#endif

73
common/Communication.h
View File

@ -244,24 +244,25 @@ inline void CommunicateMeshHalo(DoubleArray &Mesh, MPI_Comm Communicator,
{
int sendtag, recvtag;
sendtag = recvtag = 7;
PackMeshData(sendList_x, sendCount_x ,sendbuf_x, Mesh.data);
PackMeshData(sendList_X, sendCount_X ,sendbuf_X, Mesh.data);
PackMeshData(sendList_y, sendCount_y ,sendbuf_y, Mesh.data);
PackMeshData(sendList_Y, sendCount_Y ,sendbuf_Y, Mesh.data);
PackMeshData(sendList_z, sendCount_z ,sendbuf_z, Mesh.data);
PackMeshData(sendList_Z, sendCount_Z ,sendbuf_Z, Mesh.data);
PackMeshData(sendList_xy, sendCount_xy ,sendbuf_xy, Mesh.data);
PackMeshData(sendList_Xy, sendCount_Xy ,sendbuf_Xy, Mesh.data);
PackMeshData(sendList_xY, sendCount_xY ,sendbuf_xY, Mesh.data);
PackMeshData(sendList_XY, sendCount_XY ,sendbuf_XY, Mesh.data);
PackMeshData(sendList_xz, sendCount_xz ,sendbuf_xz, Mesh.data);
PackMeshData(sendList_Xz, sendCount_Xz ,sendbuf_Xz, Mesh.data);
PackMeshData(sendList_xZ, sendCount_xZ ,sendbuf_xZ, Mesh.data);
PackMeshData(sendList_XZ, sendCount_XZ ,sendbuf_XZ, Mesh.data);
PackMeshData(sendList_yz, sendCount_yz ,sendbuf_yz, Mesh.data);
PackMeshData(sendList_Yz, sendCount_Yz ,sendbuf_Yz, Mesh.data);
PackMeshData(sendList_yZ, sendCount_yZ ,sendbuf_yZ, Mesh.data);
PackMeshData(sendList_YZ, sendCount_YZ ,sendbuf_YZ, Mesh.data);
double *MeshData = Mesh.get();
PackMeshData(sendList_x, sendCount_x ,sendbuf_x, MeshData);
PackMeshData(sendList_X, sendCount_X ,sendbuf_X, MeshData);
PackMeshData(sendList_y, sendCount_y ,sendbuf_y, MeshData);
PackMeshData(sendList_Y, sendCount_Y ,sendbuf_Y, MeshData);
PackMeshData(sendList_z, sendCount_z ,sendbuf_z, MeshData);
PackMeshData(sendList_Z, sendCount_Z ,sendbuf_Z, MeshData);
PackMeshData(sendList_xy, sendCount_xy ,sendbuf_xy, MeshData);
PackMeshData(sendList_Xy, sendCount_Xy ,sendbuf_Xy, MeshData);
PackMeshData(sendList_xY, sendCount_xY ,sendbuf_xY, MeshData);
PackMeshData(sendList_XY, sendCount_XY ,sendbuf_XY, MeshData);
PackMeshData(sendList_xz, sendCount_xz ,sendbuf_xz, MeshData);
PackMeshData(sendList_Xz, sendCount_Xz ,sendbuf_Xz, MeshData);
PackMeshData(sendList_xZ, sendCount_xZ ,sendbuf_xZ, MeshData);
PackMeshData(sendList_XZ, sendCount_XZ ,sendbuf_XZ, MeshData);
PackMeshData(sendList_yz, sendCount_yz ,sendbuf_yz, MeshData);
PackMeshData(sendList_Yz, sendCount_Yz ,sendbuf_Yz, MeshData);
PackMeshData(sendList_yZ, sendCount_yZ ,sendbuf_yZ, MeshData);
PackMeshData(sendList_YZ, sendCount_YZ ,sendbuf_YZ, MeshData);
//......................................................................................
MPI_Sendrecv(sendbuf_x,sendCount_x,MPI_DOUBLE,rank_x,sendtag,
recvbuf_X,recvCount_X,MPI_DOUBLE,rank_X,recvtag,Communicator,MPI_STATUS_IGNORE);
@ -300,24 +301,24 @@ inline void CommunicateMeshHalo(DoubleArray &Mesh, MPI_Comm Communicator,
MPI_Sendrecv(sendbuf_yZ,sendCount_yZ,MPI_DOUBLE,rank_yZ,sendtag,
recvbuf_Yz,recvCount_Yz,MPI_DOUBLE,rank_Yz,recvtag,Communicator,MPI_STATUS_IGNORE);
//........................................................................................
UnpackMeshData(recvList_x, recvCount_x ,recvbuf_x, Mesh.data);
UnpackMeshData(recvList_X, recvCount_X ,recvbuf_X, Mesh.data);
UnpackMeshData(recvList_y, recvCount_y ,recvbuf_y, Mesh.data);
UnpackMeshData(recvList_Y, recvCount_Y ,recvbuf_Y, Mesh.data);
UnpackMeshData(recvList_z, recvCount_z ,recvbuf_z, Mesh.data);
UnpackMeshData(recvList_Z, recvCount_Z ,recvbuf_Z, Mesh.data);
UnpackMeshData(recvList_xy, recvCount_xy ,recvbuf_xy, Mesh.data);
UnpackMeshData(recvList_Xy, recvCount_Xy ,recvbuf_Xy, Mesh.data);
UnpackMeshData(recvList_xY, recvCount_xY ,recvbuf_xY, Mesh.data);
UnpackMeshData(recvList_XY, recvCount_XY ,recvbuf_XY, Mesh.data);
UnpackMeshData(recvList_xz, recvCount_xz ,recvbuf_xz, Mesh.data);
UnpackMeshData(recvList_Xz, recvCount_Xz ,recvbuf_Xz, Mesh.data);
UnpackMeshData(recvList_xZ, recvCount_xZ ,recvbuf_xZ, Mesh.data);
UnpackMeshData(recvList_XZ, recvCount_XZ ,recvbuf_XZ, Mesh.data);
UnpackMeshData(recvList_yz, recvCount_yz ,recvbuf_yz, Mesh.data);
UnpackMeshData(recvList_Yz, recvCount_Yz ,recvbuf_Yz, Mesh.data);
UnpackMeshData(recvList_yZ, recvCount_yZ ,recvbuf_yZ, Mesh.data);
UnpackMeshData(recvList_YZ, recvCount_YZ ,recvbuf_YZ, Mesh.data);
UnpackMeshData(recvList_x, recvCount_x ,recvbuf_x, MeshData);
UnpackMeshData(recvList_X, recvCount_X ,recvbuf_X, MeshData);
UnpackMeshData(recvList_y, recvCount_y ,recvbuf_y, MeshData);
UnpackMeshData(recvList_Y, recvCount_Y ,recvbuf_Y, MeshData);
UnpackMeshData(recvList_z, recvCount_z ,recvbuf_z, MeshData);
UnpackMeshData(recvList_Z, recvCount_Z ,recvbuf_Z, MeshData);
UnpackMeshData(recvList_xy, recvCount_xy ,recvbuf_xy, MeshData);
UnpackMeshData(recvList_Xy, recvCount_Xy ,recvbuf_Xy, MeshData);
UnpackMeshData(recvList_xY, recvCount_xY ,recvbuf_xY, MeshData);
UnpackMeshData(recvList_XY, recvCount_XY ,recvbuf_XY, MeshData);
UnpackMeshData(recvList_xz, recvCount_xz ,recvbuf_xz, MeshData);
UnpackMeshData(recvList_Xz, recvCount_Xz ,recvbuf_Xz, MeshData);
UnpackMeshData(recvList_xZ, recvCount_xZ ,recvbuf_xZ, MeshData);
UnpackMeshData(recvList_XZ, recvCount_XZ ,recvbuf_XZ, MeshData);
UnpackMeshData(recvList_yz, recvCount_yz ,recvbuf_yz, MeshData);
UnpackMeshData(recvList_Yz, recvCount_Yz ,recvbuf_Yz, MeshData);
UnpackMeshData(recvList_yZ, recvCount_yZ ,recvbuf_yZ, MeshData);
UnpackMeshData(recvList_YZ, recvCount_YZ ,recvbuf_YZ, MeshData);
}

150
common/Domain.h
View File

@ -31,8 +31,8 @@ struct Domain{
nprocx=npx; nprocy=npy; nprocz=npz;
N = Nx*Ny*Nz;
id = new char [N];
BlobLabel.New(Nx,Ny,Nz);
BlobGraph.New(18,MAX_BLOB_COUNT,MAX_BLOB_COUNT);
BlobLabel.resize(Nx,Ny,Nz);
BlobGraph.resize(18,MAX_BLOB_COUNT,MAX_BLOB_COUNT);
BoundaryCondition = BC;
}
~Domain();
@ -616,24 +616,25 @@ inline void Domain::CommunicateMeshHalo(DoubleArray &Mesh)
{
int sendtag, recvtag;
sendtag = recvtag = 7;
PackMeshData(sendList_x, sendCount_x ,sendData_x, Mesh.data);
PackMeshData(sendList_X, sendCount_X ,sendData_X, Mesh.data);
PackMeshData(sendList_y, sendCount_y ,sendData_y, Mesh.data);
PackMeshData(sendList_Y, sendCount_Y ,sendData_Y, Mesh.data);
PackMeshData(sendList_z, sendCount_z ,sendData_z, Mesh.data);
PackMeshData(sendList_Z, sendCount_Z ,sendData_Z, Mesh.data);
PackMeshData(sendList_xy, sendCount_xy ,sendData_xy, Mesh.data);
PackMeshData(sendList_Xy, sendCount_Xy ,sendData_Xy, Mesh.data);
PackMeshData(sendList_xY, sendCount_xY ,sendData_xY, Mesh.data);
PackMeshData(sendList_XY, sendCount_XY ,sendData_XY, Mesh.data);
PackMeshData(sendList_xz, sendCount_xz ,sendData_xz, Mesh.data);
PackMeshData(sendList_Xz, sendCount_Xz ,sendData_Xz, Mesh.data);
PackMeshData(sendList_xZ, sendCount_xZ ,sendData_xZ, Mesh.data);
PackMeshData(sendList_XZ, sendCount_XZ ,sendData_XZ, Mesh.data);
PackMeshData(sendList_yz, sendCount_yz ,sendData_yz, Mesh.data);
PackMeshData(sendList_Yz, sendCount_Yz ,sendData_Yz, Mesh.data);
PackMeshData(sendList_yZ, sendCount_yZ ,sendData_yZ, Mesh.data);
PackMeshData(sendList_YZ, sendCount_YZ ,sendData_YZ, Mesh.data);
double *MeshData = Mesh.get();
PackMeshData(sendList_x, sendCount_x ,sendData_x, MeshData);
PackMeshData(sendList_X, sendCount_X ,sendData_X, MeshData);
PackMeshData(sendList_y, sendCount_y ,sendData_y, MeshData);
PackMeshData(sendList_Y, sendCount_Y ,sendData_Y, MeshData);
PackMeshData(sendList_z, sendCount_z ,sendData_z, MeshData);
PackMeshData(sendList_Z, sendCount_Z ,sendData_Z, MeshData);
PackMeshData(sendList_xy, sendCount_xy ,sendData_xy, MeshData);
PackMeshData(sendList_Xy, sendCount_Xy ,sendData_Xy, MeshData);
PackMeshData(sendList_xY, sendCount_xY ,sendData_xY, MeshData);
PackMeshData(sendList_XY, sendCount_XY ,sendData_XY, MeshData);
PackMeshData(sendList_xz, sendCount_xz ,sendData_xz, MeshData);
PackMeshData(sendList_Xz, sendCount_Xz ,sendData_Xz, MeshData);
PackMeshData(sendList_xZ, sendCount_xZ ,sendData_xZ, MeshData);
PackMeshData(sendList_XZ, sendCount_XZ ,sendData_XZ, MeshData);
PackMeshData(sendList_yz, sendCount_yz ,sendData_yz, MeshData);
PackMeshData(sendList_Yz, sendCount_Yz ,sendData_Yz, MeshData);
PackMeshData(sendList_yZ, sendCount_yZ ,sendData_yZ, MeshData);
PackMeshData(sendList_YZ, sendCount_YZ ,sendData_YZ, MeshData);
//......................................................................................
MPI_Sendrecv(sendData_x,sendCount_x,MPI_DOUBLE,rank_x,sendtag,
recvData_X,recvCount_X,MPI_DOUBLE,rank_X,recvtag,Comm,MPI_STATUS_IGNORE);
@ -672,24 +673,24 @@ inline void Domain::CommunicateMeshHalo(DoubleArray &Mesh)
MPI_Sendrecv(sendData_yZ,sendCount_yZ,MPI_DOUBLE,rank_yZ,sendtag,
recvData_Yz,recvCount_Yz,MPI_DOUBLE,rank_Yz,recvtag,Comm,MPI_STATUS_IGNORE);
//........................................................................................
UnpackMeshData(recvList_x, recvCount_x ,recvData_x, Mesh.data);
UnpackMeshData(recvList_X, recvCount_X ,recvData_X, Mesh.data);
UnpackMeshData(recvList_y, recvCount_y ,recvData_y, Mesh.data);
UnpackMeshData(recvList_Y, recvCount_Y ,recvData_Y, Mesh.data);
UnpackMeshData(recvList_z, recvCount_z ,recvData_z, Mesh.data);
UnpackMeshData(recvList_Z, recvCount_Z ,recvData_Z, Mesh.data);
UnpackMeshData(recvList_xy, recvCount_xy ,recvData_xy, Mesh.data);
UnpackMeshData(recvList_Xy, recvCount_Xy ,recvData_Xy, Mesh.data);
UnpackMeshData(recvList_xY, recvCount_xY ,recvData_xY, Mesh.data);
UnpackMeshData(recvList_XY, recvCount_XY ,recvData_XY, Mesh.data);
UnpackMeshData(recvList_xz, recvCount_xz ,recvData_xz, Mesh.data);
UnpackMeshData(recvList_Xz, recvCount_Xz ,recvData_Xz, Mesh.data);
UnpackMeshData(recvList_xZ, recvCount_xZ ,recvData_xZ, Mesh.data);
UnpackMeshData(recvList_XZ, recvCount_XZ ,recvData_XZ, Mesh.data);
UnpackMeshData(recvList_yz, recvCount_yz ,recvData_yz, Mesh.data);
UnpackMeshData(recvList_Yz, recvCount_Yz ,recvData_Yz, Mesh.data);
UnpackMeshData(recvList_yZ, recvCount_yZ ,recvData_yZ, Mesh.data);
UnpackMeshData(recvList_YZ, recvCount_YZ ,recvData_YZ, Mesh.data);
UnpackMeshData(recvList_x, recvCount_x ,recvData_x, MeshData);
UnpackMeshData(recvList_X, recvCount_X ,recvData_X, MeshData);
UnpackMeshData(recvList_y, recvCount_y ,recvData_y, MeshData);
UnpackMeshData(recvList_Y, recvCount_Y ,recvData_Y, MeshData);
UnpackMeshData(recvList_z, recvCount_z ,recvData_z, MeshData);
UnpackMeshData(recvList_Z, recvCount_Z ,recvData_Z, MeshData);
UnpackMeshData(recvList_xy, recvCount_xy ,recvData_xy, MeshData);
UnpackMeshData(recvList_Xy, recvCount_Xy ,recvData_Xy, MeshData);
UnpackMeshData(recvList_xY, recvCount_xY ,recvData_xY, MeshData);
UnpackMeshData(recvList_XY, recvCount_XY ,recvData_XY, MeshData);
UnpackMeshData(recvList_xz, recvCount_xz ,recvData_xz, MeshData);
UnpackMeshData(recvList_Xz, recvCount_Xz ,recvData_Xz, MeshData);
UnpackMeshData(recvList_xZ, recvCount_xZ ,recvData_xZ, MeshData);
UnpackMeshData(recvList_XZ, recvCount_XZ ,recvData_XZ, MeshData);
UnpackMeshData(recvList_yz, recvCount_yz ,recvData_yz, MeshData);
UnpackMeshData(recvList_Yz, recvCount_Yz ,recvData_Yz, MeshData);
UnpackMeshData(recvList_yZ, recvCount_yZ ,recvData_yZ, MeshData);
UnpackMeshData(recvList_YZ, recvCount_YZ ,recvData_YZ, MeshData);
}
void Domain::BlobComm(MPI_Comm Communicator){
@ -697,24 +698,25 @@ void Domain::BlobComm(MPI_Comm Communicator){
int sendtag, recvtag;
sendtag = recvtag = 51;
//......................................................................................
PackBlobData(sendList_x, sendCount_x ,sendBuf_x, BlobLabel.data);
PackBlobData(sendList_X, sendCount_X ,sendBuf_X, BlobLabel.data);
PackBlobData(sendList_y, sendCount_y ,sendBuf_y, BlobLabel.data);
PackBlobData(sendList_Y, sendCount_Y ,sendBuf_Y, BlobLabel.data);
PackBlobData(sendList_z, sendCount_z ,sendBuf_z, BlobLabel.data);
PackBlobData(sendList_Z, sendCount_Z ,sendBuf_Z, BlobLabel.data);
PackBlobData(sendList_xy, sendCount_xy ,sendBuf_xy, BlobLabel.data);
PackBlobData(sendList_Xy, sendCount_Xy ,sendBuf_Xy, BlobLabel.data);
PackBlobData(sendList_xY, sendCount_xY ,sendBuf_xY, BlobLabel.data);
PackBlobData(sendList_XY, sendCount_XY ,sendBuf_XY, BlobLabel.data);
PackBlobData(sendList_xz, sendCount_xz ,sendBuf_xz, BlobLabel.data);
PackBlobData(sendList_Xz, sendCount_Xz ,sendBuf_Xz, BlobLabel.data);
PackBlobData(sendList_xZ, sendCount_xZ ,sendBuf_xZ, BlobLabel.data);
PackBlobData(sendList_XZ, sendCount_XZ ,sendBuf_XZ, BlobLabel.data);
PackBlobData(sendList_yz, sendCount_yz ,sendBuf_yz, BlobLabel.data);
PackBlobData(sendList_Yz, sendCount_Yz ,sendBuf_Yz, BlobLabel.data);
PackBlobData(sendList_yZ, sendCount_yZ ,sendBuf_yZ, BlobLabel.data);
PackBlobData(sendList_YZ, sendCount_YZ ,sendBuf_YZ, BlobLabel.data);
int *BlobLabelData = BlobLabel.get();
PackBlobData(sendList_x, sendCount_x ,sendBuf_x, BlobLabelData);
PackBlobData(sendList_X, sendCount_X ,sendBuf_X, BlobLabelData);
PackBlobData(sendList_y, sendCount_y ,sendBuf_y, BlobLabelData);
PackBlobData(sendList_Y, sendCount_Y ,sendBuf_Y, BlobLabelData);
PackBlobData(sendList_z, sendCount_z ,sendBuf_z, BlobLabelData);
PackBlobData(sendList_Z, sendCount_Z ,sendBuf_Z, BlobLabelData);
PackBlobData(sendList_xy, sendCount_xy ,sendBuf_xy, BlobLabelData);
PackBlobData(sendList_Xy, sendCount_Xy ,sendBuf_Xy, BlobLabelData);
PackBlobData(sendList_xY, sendCount_xY ,sendBuf_xY, BlobLabelData);
PackBlobData(sendList_XY, sendCount_XY ,sendBuf_XY, BlobLabelData);
PackBlobData(sendList_xz, sendCount_xz ,sendBuf_xz, BlobLabelData);
PackBlobData(sendList_Xz, sendCount_Xz ,sendBuf_Xz, BlobLabelData);
PackBlobData(sendList_xZ, sendCount_xZ ,sendBuf_xZ, BlobLabelData);
PackBlobData(sendList_XZ, sendCount_XZ ,sendBuf_XZ, BlobLabelData);
PackBlobData(sendList_yz, sendCount_yz ,sendBuf_yz, BlobLabelData);
PackBlobData(sendList_Yz, sendCount_Yz ,sendBuf_Yz, BlobLabelData);
PackBlobData(sendList_yZ, sendCount_yZ ,sendBuf_yZ, BlobLabelData);
PackBlobData(sendList_YZ, sendCount_YZ ,sendBuf_YZ, BlobLabelData);
//......................................................................................
MPI_Sendrecv(sendBuf_x,sendCount_x,MPI_INT,rank_x,sendtag,
recvBuf_X,recvCount_X,MPI_INT,rank_X,recvtag,Comm,MPI_STATUS_IGNORE);
@ -753,24 +755,24 @@ void Domain::BlobComm(MPI_Comm Communicator){
MPI_Sendrecv(sendBuf_yZ,sendCount_yZ,MPI_INT,rank_yZ,sendtag,
recvBuf_Yz,recvCount_Yz,MPI_INT,rank_Yz,recvtag,Comm,MPI_STATUS_IGNORE);
//........................................................................................
UnpackBlobData(recvList_x, recvCount_x ,recvBuf_x, BlobLabel.data);
UnpackBlobData(recvList_X, recvCount_X ,recvBuf_X, BlobLabel.data);
UnpackBlobData(recvList_y, recvCount_y ,recvBuf_y, BlobLabel.data);
UnpackBlobData(recvList_Y, recvCount_Y ,recvBuf_Y, BlobLabel.data);
UnpackBlobData(recvList_z, recvCount_z ,recvBuf_z, BlobLabel.data);
UnpackBlobData(recvList_Z, recvCount_Z ,recvBuf_Z, BlobLabel.data);
UnpackBlobData(recvList_xy, recvCount_xy ,recvBuf_xy, BlobLabel.data);
UnpackBlobData(recvList_Xy, recvCount_Xy ,recvBuf_Xy, BlobLabel.data);
UnpackBlobData(recvList_xY, recvCount_xY ,recvBuf_xY, BlobLabel.data);
UnpackBlobData(recvList_XY, recvCount_XY ,recvBuf_XY, BlobLabel.data);
UnpackBlobData(recvList_xz, recvCount_xz ,recvBuf_xz, BlobLabel.data);
UnpackBlobData(recvList_Xz, recvCount_Xz ,recvBuf_Xz, BlobLabel.data);
UnpackBlobData(recvList_xZ, recvCount_xZ ,recvBuf_xZ, BlobLabel.data);
UnpackBlobData(recvList_XZ, recvCount_XZ ,recvBuf_XZ, BlobLabel.data);
UnpackBlobData(recvList_yz, recvCount_yz ,recvBuf_yz, BlobLabel.data);
UnpackBlobData(recvList_Yz, recvCount_Yz ,recvBuf_Yz, BlobLabel.data);
UnpackBlobData(recvList_yZ, recvCount_yZ ,recvBuf_yZ, BlobLabel.data);
UnpackBlobData(recvList_YZ, recvCount_YZ ,recvBuf_YZ, BlobLabel.data);
UnpackBlobData(recvList_x, recvCount_x ,recvBuf_x, BlobLabelData);
UnpackBlobData(recvList_X, recvCount_X ,recvBuf_X, BlobLabelData);
UnpackBlobData(recvList_y, recvCount_y ,recvBuf_y, BlobLabelData);
UnpackBlobData(recvList_Y, recvCount_Y ,recvBuf_Y, BlobLabelData);
UnpackBlobData(recvList_z, recvCount_z ,recvBuf_z, BlobLabelData);
UnpackBlobData(recvList_Z, recvCount_Z ,recvBuf_Z, BlobLabelData);
UnpackBlobData(recvList_xy, recvCount_xy ,recvBuf_xy, BlobLabelData);
UnpackBlobData(recvList_Xy, recvCount_Xy ,recvBuf_Xy, BlobLabelData);
UnpackBlobData(recvList_xY, recvCount_xY ,recvBuf_xY, BlobLabelData);
UnpackBlobData(recvList_XY, recvCount_XY ,recvBuf_XY, BlobLabelData);
UnpackBlobData(recvList_xz, recvCount_xz ,recvBuf_xz, BlobLabelData);
UnpackBlobData(recvList_Xz, recvCount_Xz ,recvBuf_Xz, BlobLabelData);
UnpackBlobData(recvList_xZ, recvCount_xZ ,recvBuf_xZ, BlobLabelData);
UnpackBlobData(recvList_XZ, recvCount_XZ ,recvBuf_XZ, BlobLabelData);
UnpackBlobData(recvList_yz, recvCount_yz ,recvBuf_yz, BlobLabelData);
UnpackBlobData(recvList_Yz, recvCount_Yz ,recvBuf_Yz, BlobLabelData);
UnpackBlobData(recvList_yZ, recvCount_yZ ,recvBuf_yZ, BlobLabelData);
UnpackBlobData(recvList_YZ, recvCount_YZ ,recvBuf_YZ, BlobLabelData);
//......................................................................................
}

142
common/TwoPhase.h
View File

@ -201,37 +201,37 @@ public:
Volume=(Nx-2)*(Ny-2)*(Nz-2)*Dm.nprocx*Dm.nprocy*Dm.nprocz*1.0;
ncubes=(Nx-2)*(Ny-2)*(Nz-2);
cubeList.New(3,ncubes);
cubeList.resize(3,ncubes);
// Global arrays
BlobLabel.New(Nx,Ny,Nz);
SDn.New(Nx,Ny,Nz);
SDs.New(Nx,Ny,Nz);
Phase.New(Nx,Ny,Nz);
Press.New(Nx,Ny,Nz);
dPdt.New(Nx,Ny,Nz);
MeanCurvature.New(Nx,Ny,Nz);
GaussCurvature.New(Nx,Ny,Nz);
SDs_x.New(Nx,Ny,Nz); // Gradient of the signed distance
SDs_y.New(Nx,Ny,Nz);
SDs_z.New(Nx,Ny,Nz);
SDn_x.New(Nx,Ny,Nz); // Gradient of the signed distance
SDn_y.New(Nx,Ny,Nz);
SDn_z.New(Nx,Ny,Nz);
DelPhi.New(Nx,Ny,Nz);
Phase_tplus.New(Nx,Ny,Nz);
Phase_tminus.New(Nx,Ny,Nz);
Vel_x.New(Nx,Ny,Nz); // Gradient of the phase indicator field
Vel_y.New(Nx,Ny,Nz);
Vel_z.New(Nx,Ny,Nz);
BlobLabel.resize(Nx,Ny,Nz);
SDn.resize(Nx,Ny,Nz);
SDs.resize(Nx,Ny,Nz);
Phase.resize(Nx,Ny,Nz);
Press.resize(Nx,Ny,Nz);
dPdt.resize(Nx,Ny,Nz);
MeanCurvature.resize(Nx,Ny,Nz);
GaussCurvature.resize(Nx,Ny,Nz);
SDs_x.resize(Nx,Ny,Nz); // Gradient of the signed distance
SDs_y.resize(Nx,Ny,Nz);
SDs_z.resize(Nx,Ny,Nz);
SDn_x.resize(Nx,Ny,Nz); // Gradient of the signed distance
SDn_y.resize(Nx,Ny,Nz);
SDn_z.resize(Nx,Ny,Nz);
DelPhi.resize(Nx,Ny,Nz);
Phase_tplus.resize(Nx,Ny,Nz);
Phase_tminus.resize(Nx,Ny,Nz);
Vel_x.resize(Nx,Ny,Nz); // Gradient of the phase indicator field
Vel_y.resize(Nx,Ny,Nz);
Vel_z.resize(Nx,Ny,Nz);
//.........................................
// Allocate cube storage space
CubeValues.New(2,2,2);
nw_tris.New(3,20);
ns_tris.New(3,20);
ws_tris.New(3,20);
nws_seg.New(2,20);
local_sol_tris.New(3,18);
CubeValues.resize(2,2,2);
nw_tris.resize(3,20);
ns_tris.resize(3,20);
ws_tris.resize(3,20);
nws_seg.resize(2,20);
local_sol_tris.resize(3,18);
nw_pts=DTMutableList<Point>(20);
ns_pts=DTMutableList<Point>(20);
ws_pts=DTMutableList<Point>(20);
@ -240,27 +240,27 @@ public:
local_sol_pts=DTMutableList<Point>(20);
tmp=DTMutableList<Point>(20);
//.........................................
Values.New(20);
DistanceValues.New(20);
KGwns_values.New(20);
KNwns_values.New(20);
InterfaceSpeed.New(20);
NormalVector.New(60);
Values.resize(20);
DistanceValues.resize(20);
KGwns_values.resize(20);
KNwns_values.resize(20);
InterfaceSpeed.resize(20);
NormalVector.resize(60);
//.........................................
van.New(3);
vaw.New(3);
vawn.New(3);
vawns.New(3);
Gwn.New(6);
Gns.New(6);
Gws.New(6);
van_global.New(3);
vaw_global.New(3);
vawn_global.New(3);
vawns_global.New(3);
Gwn_global.New(6);
Gns_global.New(6);
Gws_global.New(6);
van.resize(3);
vaw.resize(3);
vawn.resize(3);
vawns.resize(3);
Gwn.resize(6);
Gns.resize(6);
Gws.resize(6);
van_global.resize(3);
vaw_global.resize(3);
vawn_global.resize(3);
vawns_global.resize(3);
Gwn_global.resize(6);
Gns_global.resize(6);
Gws_global.resize(6);
//.........................................
if (Dm.rank==0){
TIMELOG = fopen("timelog.tcat","a+");
@ -453,27 +453,27 @@ void TwoPhase::ComputeLocal(){
if ( Phase(i+cube[p][0],j+cube[p][1],k+cube[p][2]) > 0 ){
nwp_volume += 0.125;
// volume the excludes the interfacial region
if (DelPhi.data[n] < 1e-4){
if (DelPhi(n) < 1e-4){
vol_n += 0.125;
// pressure
pan += 0.125*Press.data[n];
pan += 0.125*Press(n);
// velocity
van(0) += 0.125*Vel_x.data[n];
van(1) += 0.125*Vel_y.data[n];
van(2) += 0.125*Vel_z.data[n];
van(0) += 0.125*Vel_x(n);
van(1) += 0.125*Vel_y(n);
van(2) += 0.125*Vel_z(n);
}
}
else{
wp_volume += 0.125;
if (DelPhi.data[n] < 1e-4){
if (DelPhi(n) < 1e-4){
// volume the excludes the interfacial region
vol_w += 0.125;
// pressure
paw += 0.125*Press.data[n];
paw += 0.125*Press(n);
// velocity
vaw(0) += 0.125*Vel_x.data[n];
vaw(1) += 0.125*Vel_y.data[n];
vaw(2) += 0.125*Vel_z.data[n];
vaw(0) += 0.125*Vel_x(n);
vaw(1) += 0.125*Vel_y(n);
vaw(2) += 0.125*Vel_z(n);
}
}
}
@ -533,14 +533,14 @@ void TwoPhase::ComputeLocalBlob(){
label=0;
nblobs_global = 0;
for (n=0; n<Nx*Ny*Nz; n++){
if (label < BlobLabel.data[n]) label = BlobLabel.data[n];
if (label < BlobLabel(n)) label = BlobLabel(n);
}
MPI_Allreduce(&label,&nblobs_global,1,MPI_INT,MPI_MAX,Dm.Comm);
if (Dm.rank==0) printf("Number of blobs is %i \n",nblobs_global);
//BlobAverages.Set(nblobs_global);
BlobAverages.New(BLOB_AVG_COUNT,nblobs_global);
BlobAverages.resize(BLOB_AVG_COUNT,nblobs_global);
BlobAverages.fill(0.0);
// Perform averaging
for (int idx=0; idx<BlobAverages.m*BlobAverages.n; idx++) BlobAverages.data[idx] = 0.0;
for (int c=0;c<ncubes;c++){
// Get cube from the list
i = cubeList(0,c);
@ -563,29 +563,29 @@ void TwoPhase::ComputeLocalBlob(){
if ( Phase(i+cube[p][0],j+cube[p][1],k+cube[p][2]) > 0 ){
BlobAverages(1,label) += 0.125;
// volume the excludes the interfacial region
if (DelPhi.data[n] < 1e-4){
if (DelPhi(n) < 1e-4){
BlobAverages(0,label) += 0.125;
// pressure
BlobAverages(2,label ) += 0.125*Press.data[n];
BlobAverages(2,label ) += 0.125*Press(n);
// velocity
BlobAverages(9,label) += 0.125*Vel_x.data[n];
BlobAverages(10,label) += 0.125*Vel_y.data[n];
BlobAverages(11,label) += 0.125*Vel_z.data[n];
BlobAverages(9,label) += 0.125*Vel_x(n);
BlobAverages(10,label) += 0.125*Vel_y(n);
BlobAverages(11,label) += 0.125*Vel_z(n);
}
}
else{
wp_volume += 0.125;
if (DelPhi.data[n] < 1e-4){
if (DelPhi(n) < 1e-4){
// volume the excludes the interfacial region
vol_w += 0.125;
// pressure
if (isnan(Press.data[n])) printf("Pressure is nan!\n");
else paw += 0.125*Press.data[n];
if (isnan(Press(n))) printf("Pressure is nan!\n");
else paw += 0.125*Press(n);
// velocity
vaw(0) += 0.125*Vel_x.data[n];
vaw(1) += 0.125*Vel_y.data[n];
vaw(2) += 0.125*Vel_z.data[n];
vaw(0) += 0.125*Vel_x(n);
vaw(1) += 0.125*Vel_y(n);
vaw(2) += 0.125*Vel_z(n);
}
}
}

8
common/pmmc.h
View File

@ -314,7 +314,7 @@ public:
DoubleArray Corners;
TriLinPoly(){
Corners.New(2,2,2);
Corners.resize(2,2,2);
}
~TriLinPoly(){
}
@ -378,9 +378,9 @@ inline int ComputeBlob(IntArray &blobs, int &nblobs, int &ncubes, IntArray &indi
// Compute the blob (F>vf|S>vs) starting from (i,j,k) - oil blob
// F>vf => oil phase S>vs => in porespace
// update the list of blobs, indicator mesh
int m = F.m; // maxima for the meshes
int n = F.n;
int o = F.o;
int m = F.size(0); // maxima for the meshes
int n = F.size(1);
int o = F.size(2);
int cubes_in_blob=0;
int nrecent = 1; // number of nodes added at most recent sweep

46
tests/BasicSimulator.cpp
View File

@ -390,10 +390,10 @@ int main(int argc, char **argv)
// sprintf(LocalRankString,"%05d",rank);
// sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
//.......................................................................
SignedDistance(SignDist.data,nspheres,cx,cy,cz,rad,Lx,Ly,Lz,Nx,Ny,Nz,
SignedDistance(SignDist.get(),nspheres,cx,cy,cz,rad,Lx,Ly,Lz,Nx,Ny,Nz,
iproc,jproc,kproc,nprocx,nprocy,nprocz);
// for (n=0; n<Nx*Ny*Nz; n++) SignDist.data[n] += (1.0); // map by a pixel to account for interface width
// for (n=0; n<Nx*Ny*Nz; n++) SignDist.get()[n] += (1.0); // map by a pixel to account for interface width
//.......................................................................
// Assign the phase ID field based on the signed distance
@ -411,11 +411,11 @@ int main(int argc, char **argv)
for ( j=1;j<Ny-1;j++){
for ( i=1;i<Nx-1;i++){
n = k*Nx*Ny+j*Nx+i;
if (SignDist.data[n] > 0.0){
if (SignDist.get()[n] > 0.0){
id[n] = 2;
}
// compute the porosity (actual interface location used)
if (SignDist.data[n] > 0.0){
if (SignDist.get()[n] > 0.0){
sum++;
}
}
@ -455,7 +455,7 @@ int main(int argc, char **argv)
for (i=0;i<Nx;i++){
n = k*Nx*Ny+j*Nx+i;
id[n] = 1;
SignDist.data[n] = max(SignDist.data[n],1.0*(2.5-k));
SignDist.get()[n] = max(SignDist.get()[n],1.0*(2.5-k));
}
}
}
@ -466,7 +466,7 @@ int main(int argc, char **argv)
for (i=0;i<Nx;i++){
n = k*Nx*Ny+j*Nx+i;
id[n] = 2;
SignDist.data[n] = max(SignDist.data[n],1.0*(k-Nz+2.5));
SignDist.get()[n] = max(SignDist.get()[n],1.0*(k-Nz+2.5));
}
}
}
@ -1036,7 +1036,7 @@ int main(int argc, char **argv)
AllocateDeviceMemory((void **) &Velocity, 3*dist_mem_size);
AllocateDeviceMemory((void **) &ColorGrad, 3*dist_mem_size);
// Copy signed distance for device initialization
CopyToDevice(dvcSignDist, SignDist.data, dist_mem_size);
CopyToDevice(dvcSignDist, SignDist.get(), dist_mem_size);
//...........................................................................
// Phase indicator (in array form as needed by PMMC algorithm)
DoubleArray Phase(Nx,Ny,Nz);
@ -1204,7 +1204,7 @@ int main(int argc, char **argv)
// sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
// WriteLocalSolidID(LocalRankFilename, id, N);
sprintf(LocalRankFilename,"%s%s","SignDist.",LocalRankString);
WriteLocalSolidDistance(LocalRankFilename, SignDist.data, N);
WriteLocalSolidDistance(LocalRankFilename, SignDist.get(), N);
//.......................................................................
if (Restart == true){
if (rank==0) printf("Reading restart file! \n");
@ -1223,7 +1223,7 @@ int main(int argc, char **argv)
InitD3Q7(ID, B_even, B_odd, &Den[N], Nx, Ny, Nz);
// Once phase has been initialized, map solid to account for 'smeared' interface
//......................................................................
for (i=0; i<N; i++) SignDist.data[i] -= (1.0); //
for (i=0; i<N; i++) SignDist.get()[i] -= (1.0); //
//.......................................................................
@ -1420,7 +1420,7 @@ int main(int argc, char **argv)
//...........................................................................
// Copy the phase indicator field for the earlier timestep
DeviceBarrier();
CopyToHost(Phase_tplus.data,Phi,N*sizeof(double));
CopyToHost(Phase_tplus.get(),Phi,N*sizeof(double));
//...........................................................................
//...........................................................................
// Copy the data for for the analysis timestep
@ -1429,11 +1429,11 @@ int main(int argc, char **argv)
//...........................................................................
DeviceBarrier();
ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
CopyToHost(Phase.data,Phi,N*sizeof(double));
CopyToHost(Press.data,Pressure,N*sizeof(double));
CopyToHost(Vel_x.data,&Velocity[0],N*sizeof(double));
CopyToHost(Vel_y.data,&Velocity[N],N*sizeof(double));
CopyToHost(Vel_z.data,&Velocity[2*N],N*sizeof(double));
CopyToHost(Phase.get(),Phi,N*sizeof(double));
CopyToHost(Press.get(),Pressure,N*sizeof(double));
CopyToHost(Vel_x.get(),&Velocity[0],N*sizeof(double));
CopyToHost(Vel_y.get(),&Velocity[N],N*sizeof(double));
CopyToHost(Vel_z.get(),&Velocity[2*N],N*sizeof(double));
MPI_Barrier(MPI_COMM_WORLD);
//...........................................................................
@ -1881,7 +1881,7 @@ int main(int argc, char **argv)
//...........................................................................
// Copy the phase indicator field for the earlier timestep
DeviceBarrier();
CopyToHost(Phase_tplus.data,Phi,N*sizeof(double));
CopyToHost(Phase_tplus.get(),Phi,N*sizeof(double));
//...........................................................................
}
if (timestep%1000 == 0){
@ -1892,18 +1892,18 @@ int main(int argc, char **argv)
//...........................................................................
DeviceBarrier();
ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
CopyToHost(Phase.data,Phi,N*sizeof(double));
CopyToHost(Press.data,Pressure,N*sizeof(double));
CopyToHost(Vel_x.data,&Velocity[0],N*sizeof(double));
CopyToHost(Vel_y.data,&Velocity[N],N*sizeof(double));
CopyToHost(Vel_z.data,&Velocity[2*N],N*sizeof(double));
CopyToHost(Phase.get(),Phi,N*sizeof(double));
CopyToHost(Press.get(),Pressure,N*sizeof(double));
CopyToHost(Vel_x.get(),&Velocity[0],N*sizeof(double));
CopyToHost(Vel_y.get(),&Velocity[N],N*sizeof(double));
CopyToHost(Vel_z.get(),&Velocity[2*N],N*sizeof(double));
MPI_Barrier(MPI_COMM_WORLD);
}
if (timestep%1000 == 5){
//...........................................................................
// Copy the phase indicator field for the later timestep
DeviceBarrier();
CopyToHost(Phase_tminus.data,Phi,N*sizeof(double));
CopyToHost(Phase_tminus.get(),Phi,N*sizeof(double));
//...........................................................................
// Calculate the time derivative of the phase indicator field
for (n=0; n<N; n++) dPdt(n) = 0.1*(Phase_tplus(n) - Phase_tminus(n));
@ -2413,7 +2413,7 @@ int main(int argc, char **argv)
sprintf(LocalRankFilename,"%s/%s%s",tmpstr,"dPdt.",LocalRankString);
SPEED = fopen(LocalRankFilename,"wb");
fwrite(dPdt.data,8,N,SPEED);
fwrite(dPdt.get(),8,N,SPEED);
fclose(SPEED);
}

31
tests/BlobAnalysis.cpp
View File

@ -437,9 +437,8 @@ int main(int argc, char **argv)
}
*/
FILE *PHASE;
PHASE = fopen("Phase.dat","wb");
fwrite(Phase.data,8,Nx*Ny*Nz,PHASE);
FILE *PHASE = fopen("Phase.dat","wb");
fwrite(Phase.get(),8,Nx*Ny*Nz,PHASE);
fclose(PHASE);
// Initialize the local blob ID
@ -590,9 +589,9 @@ int main(int argc, char **argv)
// Otherwise, this point has already been assigned - ignore
// Make sure list blob_nodes is large enough
if ( nblobs > b.Length-1){
if ( nblobs > (int)b.length()-1){
printf("Increasing size of blob list \n");
b = IncreaseSize(b,b.Length);
b.resize(2*b.length());
}
}
}
@ -633,7 +632,7 @@ int main(int argc, char **argv)
DoubleArray BlobAverages(NUM_AVERAGES,nblobs);
// Map the signed distance for the analysis
for (i=0; i<Nx*Ny*Nz; i++) SignDist.data[i] -= (1.0);
for (i=0; i<Nx*Ny*Nz; i++) SignDist(i) -= (1.0);
// Compute the porosity
porosity=0.0;
@ -731,11 +730,11 @@ int main(int argc, char **argv)
// volume the excludes the interfacial region
vol_n += 0.125;
// pressure
pan += 0.125*Press.data[n];
pan += 0.125*Press(n);
// velocity
van(0) += 0.125*Vel_x.data[n];
van(1) += 0.125*Vel_y.data[n];
van(2) += 0.125*Vel_z.data[n];
van(0) += 0.125*Vel_x(n);
van(1) += 0.125*Vel_y(n);
van(2) += 0.125*Vel_z(n);
}
// volume averages over the wetting phase
@ -743,11 +742,11 @@ int main(int argc, char **argv)
// volume the excludes the interfacial region
vol_w += 0.125;
// pressure
paw += 0.125*Press.data[n];
paw += 0.125*Press(n);
// velocity
vaw(0) += 0.125*Vel_x.data[n];
vaw(1) += 0.125*Vel_y.data[n];
vaw(2) += 0.125*Vel_z.data[n];
vaw(0) += 0.125*Vel_x(n);
vaw(1) += 0.125*Vel_y(n);
vaw(2) += 0.125*Vel_z(n);
}
}
}
@ -1076,12 +1075,12 @@ int main(int argc, char **argv)
FILE *BLOBS;
BLOBS = fopen("Blobs.dat","wb");
fwrite(LocalBlobID.data,4,Nx*Ny*Nz,BLOBS);
fwrite(LocalBlobID.get(),4,Nx*Ny*Nz,BLOBS);
fclose(BLOBS);
FILE *DISTANCE;
DISTANCE = fopen("SignDist.dat","wb");
fwrite(SignDist.data,8,Nx*Ny*Nz,DISTANCE);
fwrite(SignDist.get(),8,Nx*Ny*Nz,DISTANCE);
fclose(DISTANCE);
}

6
tests/BlobIdentify.cpp
View File

@ -343,9 +343,9 @@ int main(int argc, char **argv)
// Otherwise, this point has already been assigned - ignore
// Make sure list blob_nodes is large enough
if ( nblobs > b.Length-1){
if ( nblobs > (int)b.length()-1){
printf("Increasing size of blob list \n");
b = IncreaseSize(b,b.Length);
b.resize(2*b.length());
}
}
}
@ -433,7 +433,7 @@ int main(int argc, char **argv)
FILE *BLOBS;
BLOBS = fopen("Blobs.dat","wb");
fwrite(GlobalBlobID.data,4,Nx*Ny*Nz,BLOBS);
fwrite(GlobalBlobID.get(),4,Nx*Ny*Nz,BLOBS);
fclose(BLOBS);
}

38
tests/TestBubble.cpp
View File

@ -1193,21 +1193,21 @@ int main(int argc, char **argv)
//...........................................................................
InitD3Q19(ID, f_even, f_odd, Nx, Ny, Nz);
//......................................................................
// InitDenColorDistance(ID, Copy, Phi, SDs.data, das, dbs, beta, xIntPos, Nx, Ny, Nz, S);
InitDenColorDistance(ID, Den, Phi, SDs.data, das, dbs, beta, xIntPos, Nx, Ny, Nz);
// InitDenColorDistance(ID, Copy, Phi, SDs.get(), das, dbs, beta, xIntPos, Nx, Ny, Nz, S);
InitDenColorDistance(ID, Den, Phi, SDs.get(), das, dbs, beta, xIntPos, Nx, Ny, Nz);
InitD3Q7(ID, A_even, A_odd, &Den[0], Nx, Ny, Nz);
InitD3Q7(ID, B_even, B_odd, &Den[N], Nx, Ny, Nz);
//......................................................................
// Once phase has been initialized, map solid to account for 'smeared' interface
//......................................................................
for (i=0; i<N; i++) SDs.data[i] -= (1.0); //
for (i=0; i<N; i++) SDs(i) -= (1.0); //
//......................................................................
//.......................................................................
sprintf(LocalRankString,"%05d",rank);
sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
WriteLocalSolidID(LocalRankFilename, id, N);
sprintf(LocalRankFilename,"%s%s","SDs.",LocalRankString);
WriteLocalSolidDistance(LocalRankFilename, SDs.data, N);
WriteLocalSolidDistance(LocalRankFilename, SDs.get(), N);
//.......................................................................
if (Restart == true){
if (rank==0) printf("Reading restart file! \n");
@ -1340,7 +1340,7 @@ int main(int argc, char **argv)
//...........................................................................
// Copy the phase indicator field for the earlier timestep
DeviceBarrier();
CopyToHost(Phase_tplus.data,Phi,N*sizeof(double));
CopyToHost(Phase_tplus.get(),Phi,N*sizeof(double));
//...........................................................................
//...........................................................................
// Copy the data for for the analysis timestep
@ -1349,8 +1349,8 @@ int main(int argc, char **argv)
//...........................................................................
DeviceBarrier();
ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
CopyToHost(Phase.data,Phi,N*sizeof(double));
CopyToHost(Press.data,Pressure,N*sizeof(double));
CopyToHost(Phase.get(),Phi,N*sizeof(double));
CopyToHost(Press.get(),Pressure,N*sizeof(double));
MPI_Barrier(MPI_COMM_WORLD);
//...........................................................................
@ -1769,15 +1769,15 @@ int main(int argc, char **argv)
// Copy the phase indicator field for the later timestep
DeviceBarrier();
ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
CopyToHost(Phase_tminus.data,Phi,N*sizeof(double));
CopyToHost(Phase_tplus.data,Phi,N*sizeof(double));
CopyToHost(Phase.data,Phi,N*sizeof(double));
CopyToHost(Press.data,Pressure,N*sizeof(double));
CopyToHost(Phase_tminus.get(),Phi,N*sizeof(double));
CopyToHost(Phase_tplus.get(),Phi,N*sizeof(double));
CopyToHost(Phase.get(),Phi,N*sizeof(double));
CopyToHost(Press.get(),Pressure,N*sizeof(double));
double temp=0.5/beta;
for (n=0; n<N; n++){
double value = Phase.data[n];
SDn.data[n] = temp*log((1.0+value)/(1.0-value));
double value = Phase.get()[n];
SDn(n) = temp*log((1.0+value)/(1.0-value));
}
//...........................................................................
@ -1972,7 +1972,7 @@ int main(int argc, char **argv)
// 1-D index for this cube corner
n = i+cube[p][0] + (j+cube[p][1])*Nx + (k+cube[p][2])*Nx*Ny;
// compute the norm of the gradient of the phase indicator field
delphi = sqrt(Phase_x.data[n]*Phase_x.data[n]+Phase_y.data[n]*Phase_y.data[n]+Phase_z.data[n]*Phase_z.data[n]);
delphi = sqrt(Phase_x(n)*Phase_x(n)+Phase_y(n)*Phase_y(n)+Phase_z(n)*Phase_z(n));
// Compute the non-wetting phase volume contribution
if ( Phase(i+cube[p][0],j+cube[p][1],k+cube[p][2]) > 0 ){
nwp_volume += 0.125;
@ -1980,14 +1980,14 @@ int main(int argc, char **argv)
if (delphi < 1e-4){
vol_n += 0.125;
// pressure
pan += 0.125*Press.data[n];
pan += 0.125*Press(n);
}
}
else if (delphi < 1e-4){
// volume the excludes the interfacial region
vol_w += 0.125;
// pressure
paw += 0.125*Press.data[n];
paw += 0.125*Press(n);
}
}
}
@ -2212,12 +2212,12 @@ int main(int argc, char **argv)
//************************************************************************/
sprintf(LocalRankFilename,"%s%s","Phase.",LocalRankString);
// printf("Local File Name = %s \n",LocalRankFilename);
// CopyToHost(Phase.data,Phi,N*sizeof(double));
// CopyToHost(Phase.get(),Phi,N*sizeof(double));
FILE *PHASE;
PHASE = fopen(LocalRankFilename,"wb");
fwrite(Press.data,8,N,PHASE);
// fwrite(MeanCurvature.data,8,N,PHASE);
fwrite(Press.get(),8,N,PHASE);
// fwrite(MeanCurvature.get(),8,N,PHASE);
fclose(PHASE);
/* double *DensityValues;

4
tests/TestSphereCurvature.cpp
View File

@ -80,7 +80,7 @@ int main(int argc, char **argv)
DoubleArray CubeValues(2,2,2);
// Compute the signed distance function
SignedDistance(Phase.data,nspheres,cx,cy,cz,rad,Lx,Ly,Lz,Nx,Ny,Nz,0,0,0,1,1,1);
SignedDistance(Phase.get(),nspheres,cx,cy,cz,rad,Lx,Ly,Lz,Nx,Ny,Nz,0,0,0,1,1,1);
for (k=0; k<Nz; k++){
for (j=0; j<Ny; j++){
@ -90,7 +90,7 @@ int main(int argc, char **argv)
}
}
}
SignedDistance(SignDist.data,0,cx,cy,cz,rad,Lx,Ly,Lz,Nx,Ny,Nz,0,0,0,1,1,1);
SignedDistance(SignDist.get(),0,cx,cy,cz,rad,Lx,Ly,Lz,Nx,Ny,Nz,0,0,0,1,1,1);
pmmc_MeshCurvature(Phase, MeanCurvature, GaussCurvature, Nx, Ny, Nz);

2
tests/TestTwoPhase.cpp
View File

@ -110,7 +110,7 @@ int main(int argc, char **argv)
if (rank==0){
FILE *PHASE;
PHASE = fopen("Phase.00000","wb");
fwrite(Averages.MeanCurvature.data,8,Nx*Ny*Nz,PHASE);
fwrite(Averages.MeanCurvature.get(),8,Nx*Ny*Nz,PHASE);
fclose(PHASE);
}
// ****************************************************

20
tests/TestWriter.cpp
View File

@ -55,19 +55,19 @@ int main(int argc, char **argv)
};
// Create the meshes
shared_ptr<IO::PointList> set1( new IO::PointList(N_points) );
std::shared_ptr<IO::PointList> set1( new IO::PointList(N_points) );
for (int i=0; i<N_points; i++) {
set1->points[i].x = x[i];
set1->points[i].y = y[i];
set1->points[i].z = z[i];
}
shared_ptr<IO::TriMesh> trimesh( new IO::TriMesh(N_tri,set1) );
std::shared_ptr<IO::TriMesh> trimesh( new IO::TriMesh(N_tri,set1) );
for (int i=0; i<N_tri; i++) {
trimesh->A[i] = tri[i][0];
trimesh->B[i] = tri[i][1];
trimesh->C[i] = tri[i][2];
}
shared_ptr<IO::TriList> trilist( new IO::TriList(*trimesh) );
std::shared_ptr<IO::TriList> trilist( new IO::TriList(*trimesh) );
for (int i=0; i<N_tri; i++) {
Point A(x[tri[i][0]],y[tri[i][0]],z[tri[i][0]]);
Point B(x[tri[i][1]],y[tri[i][1]],z[tri[i][1]]);
@ -80,8 +80,8 @@ int main(int argc, char **argv)
}
// Create the variables
shared_ptr<IO::Variable> dist_set1( new IO::Variable() );
shared_ptr<IO::Variable> dist_list( new IO::Variable() );
std::shared_ptr<IO::Variable> dist_set1( new IO::Variable() );
std::shared_ptr<IO::Variable> dist_list( new IO::Variable() );
dist_set1->dim = 1;
dist_list->dim = 1;
dist_set1->name = "Distance";
@ -143,7 +143,7 @@ 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++) {
shared_ptr<IO::Mesh> mesh = IO::getMesh(".",timesteps[i],list[j],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());
pass = false;
@ -151,7 +151,7 @@ int main(int argc, char **argv)
}
if ( meshData[j].meshName=="pointmesh" ) {
// Check the pointmesh
shared_ptr<IO::PointList> pmesh = IO::getPointList(mesh);
std::shared_ptr<IO::PointList> pmesh = IO::getPointList(mesh);
if ( pmesh.get()==NULL ) {
pass = false;
break;
@ -163,8 +163,8 @@ int main(int argc, char **argv)
}
if ( meshData[j].meshName=="trimesh" || meshData[j].meshName=="trilist" ) {
// Check the trimesh/trilist
shared_ptr<IO::TriMesh> mesh1 = IO::getTriMesh(mesh);
shared_ptr<IO::TriList> mesh2 = IO::getTriList(mesh);
std::shared_ptr<IO::TriMesh> mesh1 = IO::getTriMesh(mesh);
std::shared_ptr<IO::TriList> mesh2 = IO::getTriList(mesh);
if ( mesh1.get()==NULL || mesh2.get()==NULL ) {
pass = false;
break;
@ -213,7 +213,7 @@ int main(int argc, char **argv)
}
for (size_t v=0; v<list[i].variables.size(); v++) {
for (size_t k=0; k<list[i].domains.size(); k++) {
shared_ptr<const IO::Variable> variable =
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];
const IO::Variable& var2 = *variable;

58
tests/lb2_Color_wia_mpi.cpp
View File

@ -400,10 +400,10 @@ int main(int argc, char **argv)
// sprintf(LocalRankString,"%05d",rank);
// sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
//.......................................................................
SignedDistance(Averages.SDs.data,nspheres,cx,cy,cz,rad,Lx,Ly,Lz,Nx,Ny,Nz,
SignedDistance(Averages.SDs.get(),nspheres,cx,cy,cz,rad,Lx,Ly,Lz,Nx,Ny,Nz,
iproc,jproc,kproc,nprocx,nprocy,nprocz);
// for (n=0; n<Nx*Ny*Nz; n++) SDs.data[n] += (1.0); // map by a pixel to account for interface width
// for (n=0; n<Nx*Ny*Nz; n++) SDs(n) += (1.0); // map by a pixel to account for interface width
//.......................................................................
// Assign the phase ID field based on the signed distance
@ -422,11 +422,11 @@ int main(int argc, char **argv)
for ( j=1;j<Ny-1;j++){
for ( i=1;i<Nx-1;i++){
n = k*Nx*Ny+j*Nx+i;
if (Averages.SDs.data[n] > 0.0){
if (Averages.SDs(n) > 0.0){
id[n] = 2;
}
// compute the porosity (actual interface location used)
if (Averages.SDs.data[n] > 0.0){
if (Averages.SDs(n) > 0.0){
sum++;
}
}
@ -465,7 +465,7 @@ int main(int argc, char **argv)
for (i=0;i<Nx;i++){
n = k*Nx*Ny+j*Nx+i;
id[n] = 1;
Averages.SDs.data[n] = max(Averages.SDs.data[n],1.0*(2.5-k));
Averages.SDs(n) = max(Averages.SDs(n),1.0*(2.5-k));
}
}
}
@ -476,7 +476,7 @@ int main(int argc, char **argv)
for (i=0;i<Nx;i++){
n = k*Nx*Ny+j*Nx+i;
id[n] = 2;
Averages.SDs.data[n] = max(Averages.SDs.data[n],1.0*(k-Nz+2.5));
Averages.SDs(n) = max(Averages.SDs(n),1.0*(k-Nz+2.5));
}
}
}
@ -1051,7 +1051,7 @@ int main(int argc, char **argv)
AllocateDeviceMemory((void **) &Velocity, 3*dist_mem_size);
AllocateDeviceMemory((void **) &ColorGrad, 3*dist_mem_size);
// Copy signed distance for device initialization
CopyToDevice(dvcSDs, Averages.SDs.data, dist_mem_size);
CopyToDevice(dvcSDs, Averages.SDs.get(), dist_mem_size);
//...........................................................................
//copies of data needed to perform checkpointing from cpu
@ -1093,13 +1093,13 @@ int main(int argc, char **argv)
// sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
// WriteLocalSolidID(LocalRankFilename, id, N);
sprintf(LocalRankFilename,"%s%s","SDs.",LocalRankString);
WriteLocalSolidDistance(LocalRankFilename, Averages.SDs.data, N);
WriteLocalSolidDistance(LocalRankFilename, Averages.SDs.get(), N);
//......................................................................
InitD3Q7(ID, A_even, A_odd, &Den[0], Nx, Ny, Nz);
InitD3Q7(ID, B_even, B_odd, &Den[N], Nx, Ny, Nz);
// Once phase has been initialized, map solid to account for 'smeared' interface
//......................................................................
for (i=0; i<N; i++) Averages.SDs.data[i] -= (1.0); //
for (i=0; i<N; i++) Averages.SDs(i) -= (1.0); //
//.......................................................................
// Finalize setup for averaging domain
Averages.SetupCubes(Dm);
@ -1226,7 +1226,7 @@ int main(int argc, char **argv)
//...........................................................................
// Copy the phase indicator field for the earlier timestep
DeviceBarrier();
CopyToHost(Averages.Phase_tplus.data,Phi,N*sizeof(double));
CopyToHost(Averages.Phase_tplus.get(),Phi,N*sizeof(double));
//...........................................................................
//...........................................................................
// Copy the data for for the analysis timestep
@ -1235,11 +1235,11 @@ int main(int argc, char **argv)
//...........................................................................
DeviceBarrier();
ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
CopyToHost(Averages.Phase.data,Phi,N*sizeof(double));
CopyToHost(Averages.Press.data,Pressure,N*sizeof(double));
CopyToHost(Averages.Vel_x.data,&Velocity[0],N*sizeof(double));
CopyToHost(Averages.Vel_y.data,&Velocity[N],N*sizeof(double));
CopyToHost(Averages.Vel_z.data,&Velocity[2*N],N*sizeof(double));
CopyToHost(Averages.Phase.get(),Phi,N*sizeof(double));
CopyToHost(Averages.Press.get(),Pressure,N*sizeof(double));
CopyToHost(Averages.Vel_x.get(),&Velocity[0],N*sizeof(double));
CopyToHost(Averages.Vel_y.get(),&Velocity[N],N*sizeof(double));
CopyToHost(Averages.Vel_z.get(),&Velocity[2*N],N*sizeof(double));
MPI_Barrier(MPI_COMM_WORLD);
//...........................................................................
@ -1659,8 +1659,8 @@ int main(int argc, char **argv)
//...........................................................................
// Copy the phase indicator field for the earlier timestep
DeviceBarrier();
CopyToHost(Averages.Phase.data,Phi,N*sizeof(double));
Averages.ColorToSignedDistance(beta,Averages.Phase.data,Averages.Phase_tplus.data);
CopyToHost(Averages.Phase.get(),Phi,N*sizeof(double));
Averages.ColorToSignedDistance(beta,Averages.Phase.get(),Averages.Phase_tplus.get());
//...........................................................................
}
if (timestep%1000 == 0){
@ -1671,23 +1671,23 @@ int main(int argc, char **argv)
//...........................................................................
DeviceBarrier();
ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
CopyToHost(Averages.Phase.data,Phi,N*sizeof(double));
CopyToHost(Averages.Press.data,Pressure,N*sizeof(double));
CopyToHost(Averages.Vel_x.data,&Velocity[0],N*sizeof(double));
CopyToHost(Averages.Vel_y.data,&Velocity[N],N*sizeof(double));
CopyToHost(Averages.Vel_z.data,&Velocity[2*N],N*sizeof(double));
CopyToHost(Averages.Phase.get(),Phi,N*sizeof(double));
CopyToHost(Averages.Press.get(),Pressure,N*sizeof(double));
CopyToHost(Averages.Vel_x.get(),&Velocity[0],N*sizeof(double));
CopyToHost(Averages.Vel_y.get(),&Velocity[N],N*sizeof(double));
CopyToHost(Averages.Vel_z.get(),&Velocity[2*N],N*sizeof(double));
MPI_Barrier(MPI_COMM_WORLD);
}
if (timestep%1000 == 5){
//...........................................................................
// Copy the phase indicator field for the later timestep
DeviceBarrier();
CopyToHost(Averages.Phase_tminus.data,Phi,N*sizeof(double));
Averages.ColorToSignedDistance(beta,Averages.Phase_tminus.data,Averages.Phase_tminus.data);
CopyToHost(Averages.Phase_tminus.get(),Phi,N*sizeof(double));
Averages.ColorToSignedDistance(beta,Averages.Phase_tminus.get(),Averages.Phase_tminus.get());
//....................................................................
Averages.Initialize();
Averages.ComputeDelPhi();
Averages.ColorToSignedDistance(beta,Averages.Phase.data,Averages.SDn.data);
Averages.ColorToSignedDistance(beta,Averages.Phase.get(),Averages.SDn.get());
Averages.UpdateMeshValues();
Averages.ComputeLocal();
Averages.Reduce();
@ -1956,20 +1956,20 @@ int main(int argc, char **argv)
sprintf(LocalRankFilename,"%s%s","Phase.",LocalRankString);
FILE *PHASE;
PHASE = fopen(LocalRankFilename,"wb");
fwrite(Averages.Phase.data,8,N,PHASE);
fwrite(Averages.Phase.get(),8,N,PHASE);
fclose(PHASE);
sprintf(LocalRankFilename,"%s%s","Pressure.",LocalRankString);
FILE *PRESS;
PRESS = fopen(LocalRankFilename,"wb");
fwrite(Averages.Press.data,8,N,PRESS);
fwrite(Averages.Press.get(),8,N,PRESS);
fclose(PRESS);
/* sprintf(LocalRankFilename,"%s%s","dPdt.",LocalRankString);
FILE *SPEED;
SPEED = fopen(LocalRankFilename,"wb");
fwrite(dPdt.data,8,N,SPEED);
fwrite(dPdt.get(),8,N,SPEED);
fclose(SPEED);
sprintf(LocalRankFilename,"%s%s","DenA.",LocalRankString);
@ -1994,7 +1994,7 @@ int main(int argc, char **argv)
}
}
}
fwrite(Phase.data,8,N,GRAD);
fwrite(Phase.get(),8,N,GRAD);
fclose(GRAD);
*/
/* double *DensityValues;

42
tests/lbpm_BlobAnalysis.cpp
View File

@ -120,16 +120,16 @@ int main(int argc, char **argv)
if (rank == 0) cout << "Reading in domain from signed distance function..." << endl;
//.......................................................................
sprintf(LocalRankFilename,"%s%s","SignDist.",LocalRankString);
ReadBinaryFile(LocalRankFilename, Averages.SDs.data, N);
ReadBinaryFile(LocalRankFilename, Averages.SDs.get(), N);
MPI_Barrier(MPI_COMM_WORLD);
// sprintf(LocalRankFilename,"%s%s","Pressure.",LocalRankString);
//ReadBinaryFile(LocalRankFilename, Averages.Press.data, N);
//ReadBinaryFile(LocalRankFilename, Averages.Press.get(), N);
//MPI_Barrier(MPI_COMM_WORLD);
if (rank == 0) cout << "Domain set." << endl;
//.......................................................................
sprintf(LocalRankFilename,"%s%s","BlobLabel.",LocalRankString);
ReadBlobFile(LocalRankFilename, Averages.BlobLabel.data, N);
ReadBlobFile(LocalRankFilename, Averages.BlobLabel.get(), N);
MPI_Barrier(MPI_COMM_WORLD);
if (rank == 0) cout << "BlobLabel set." << endl;
@ -181,14 +181,14 @@ int main(int argc, char **argv)
vx = f1-f2+f7-f8+f9-f10+f11-f12+f13-f14;
vy = f3-f4+f7-f8-f9+f10+f15-f16+f17-f18;
vz = f5-f6+f11-f12-f13+f14+f15-f16-f17+f18;
Averages.Phase.data[n]=(da-db)/(da+db);
Averages.Phase_tplus.data[n]=(da-db)/(da+db);
Averages.Phase_tminus.data[n]=(da-db)/(da+db);
Averages.Press.data[n]=press;
Averages.Vel_x.data[n]=vx;
Averages.Vel_y.data[n]=vy;
Averages.Vel_z.data[n]=vz;
if (Averages.SDs.data[n] > 0.0){
Averages.Phase(n)=(da-db)/(da+db);
Averages.Phase_tplus(n)=(da-db)/(da+db);
Averages.Phase_tminus(n)=(da-db)/(da+db);
Averages.Press(n)=press;
Averages.Vel_x(n)=vx;
Averages.Vel_y(n)=vy;
Averages.Vel_z(n)=vz;
if (Averages.SDs(n) > 0.0){
Dm.id[n]=1;
sum += 1.0;
}
@ -201,7 +201,7 @@ int main(int argc, char **argv)
porosity = sum_global/Dm.Volume;
if (rank==0) printf("Porosity = %f \n",porosity);
Dm.CommInit(MPI_COMM_WORLD);
for (int i=0; i<N; i++) Averages.SDs.data[i] -= 1.0; // map the distance
for (int i=0; i<N; i++) Averages.SDs(i) -= 1.0; // map the distance
double beta = 0.95;
@ -209,16 +209,16 @@ int main(int argc, char **argv)
Averages.UpdateSolid();
Averages.Initialize();
Averages.ComputeDelPhi();
Averages.ColorToSignedDistance(beta,Averages.Phase.data,Averages.SDn.data);
Averages.ColorToSignedDistance(beta,Averages.Phase.get(),Averages.SDn.get());
Averages.UpdateMeshValues();
Averages.ComputeLocalBlob();
Averages.Reduce();
int b=0;
// Blobs.Set(Averages.BlobAverages.NBLOBS);
int dimx = Averages.BlobAverages.m;
int dimy = Averages.BlobAverages.n;
int TotalBlobInfoSize=Averages.BlobAverages.m*Averages.BlobAverages.n;
int dimx = Averages.BlobAverages.size(0);
int dimy = Averages.BlobAverages.size(1);
int TotalBlobInfoSize=dimx*dimy;
FILE *BLOBLOG;
if (rank==0){
@ -227,11 +227,11 @@ int main(int argc, char **argv)
}
// BlobContainer Blobs;
DoubleArray RecvBuffer(dimx);
// MPI_Allreduce(&Averages.BlobAverages.data,&Blobs.data,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
// MPI_Allreduce(&Averages.BlobAverages.get(),&Blobs.get(),1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Barrier(MPI_COMM_WORLD);
if (rank==0) printf("All ranks passed gate \n");
for (int b=0; b<Averages.BlobAverages.n; b++){
for (int b=0; b<(int)Averages.BlobAverages.size(1); b++){
MPI_Allreduce(&Averages.BlobAverages(0,b),&RecvBuffer(0),dimx,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
for (int idx=0; idx<dimx-1; idx++) Averages.BlobAverages(idx,b)=RecvBuffer(idx);
@ -271,7 +271,7 @@ int main(int argc, char **argv)
if (rank==0){
// printf("Reduced blob %i \n",b);
fprintf(BLOBLOG,"%.5g %.5g %.5g\n",Averages.vol_w_global,Averages.paw_global,Averages.aws_global);
for (int b=0; b<Averages.BlobAverages.n; b++){
for (int b=0; b<(int)Averages.BlobAverages.size(1); b++){
if (Averages.BlobAverages(0,b) > 0.0){
double Vn,pn,awn,ans,Jwn,Kwn,lwns,cwns;
Vn = Averages.BlobAverages(1,b);
@ -318,7 +318,7 @@ inline void WriteBlobStates(TwoPhase TCAT, double D, double porosity){
pw = TCAT.paw_global;
aws = TCAT.aws;
// Compute the averages over the entire non-wetting phsae
for (a=0; a<TCAT.BlobAverages.n; a++){
for (a=0; a<(int)TCAT.BlobAverages.size(1); a++){
vol_n += TCAT.BlobAverages(0,a);
pan += TCAT.BlobAverages(2,a)*TCAT.BlobAverages(0,a);
awn += TCAT.BlobAverages(3,a);
@ -334,7 +334,7 @@ inline void WriteBlobStates(TwoPhase TCAT, double D, double porosity){
// Compute the pore voume (sum of wetting an non-wetting phase volumes)
PoreVolume=TCAT.wp_volume_global + nwp_volume;
// Subtract off portions of non-wetting phase in order of size
for (a=TCAT.BlobAverages.n-1; a>0; a--){
for (a=TCAT.BlobAverages.size(1)-1; a>0; a--){
// Subtract the features one-by-one
vol_n -= TCAT.BlobAverages(0,a);
pan -= TCAT.BlobAverages(2,a)*TCAT.BlobAverages(0,a);

56
tests/lbpm_color_simulator.cpp
View File

@ -326,7 +326,7 @@ int main(int argc, char **argv)
// sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
// WriteLocalSolidID(LocalRankFilename, id, N);
sprintf(LocalRankFilename,"%s%s","SignDist.",LocalRankString);
ReadBinaryFile(LocalRankFilename, Averages.SDs.data, N);
ReadBinaryFile(LocalRankFilename, Averages.SDs.get(), N);
MPI_Barrier(MPI_COMM_WORLD);
if (rank == 0) cout << "Domain set." << endl;
@ -347,11 +347,11 @@ int main(int argc, char **argv)
for ( j=1;j<Ny-1;j++){
for ( i=1;i<Nx-1;i++){
n = k*Nx*Ny+j*Nx+i;
if (Averages.SDs.data[n] > 0.0){
if (Averages.SDs(n) > 0.0){
id[n] = 2;
}
// compute the porosity (actual interface location used)
if (Averages.SDs.data[n] > 0.0){
if (Averages.SDs(n) > 0.0){
sum++;
}
}
@ -393,7 +393,7 @@ int main(int argc, char **argv)
for (i=0;i<Nx;i++){
n = k*Nx*Ny+j*Nx+i;
id[n] = 1;
Averages.SDs.data[n] = max(Averages.SDs.data[n],1.0*(2.5-k));
Averages.SDs(n) = max(Averages.SDs(n),1.0*(2.5-k));
}
}
}
@ -404,7 +404,7 @@ int main(int argc, char **argv)
for (i=0;i<Nx;i++){
n = k*Nx*Ny+j*Nx+i;
id[n] = 2;
Averages.SDs.data[n] = max(Averages.SDs.data[n],1.0*(k-Nz+2.5));
Averages.SDs(n) = max(Averages.SDs(n),1.0*(k-Nz+2.5));
}
}
}
@ -974,7 +974,7 @@ int main(int argc, char **argv)
// Copy signed distance for device initialization
CopyToDevice(dvcSignDist, Averages.SDs.data, dist_mem_size);
CopyToDevice(dvcSignDist, Averages.SDs.get(), dist_mem_size);
//...........................................................................
int logcount = 0; // number of surface write-outs
@ -1007,7 +1007,7 @@ int main(int argc, char **argv)
InitD3Q7(ID, B_even, B_odd, &Den[N], Nx, Ny, Nz);
//.......................................................................
// Once phase has been initialized, map solid to account for 'smeared' interface
for (i=0; i<N; i++) Averages.SDs.data[i] -= (1.0); //
for (i=0; i<N; i++) Averages.SDs(i) -= (1.0); //
//.......................................................................
// Finalize setup for averaging domain
Averages.SetupCubes(Dm);
@ -1127,7 +1127,7 @@ int main(int argc, char **argv)
//...........................................................................
// Copy the phase indicator field for the earlier timestep
DeviceBarrier();
CopyToHost(Averages.Phase_tplus.data,Phi,N*sizeof(double));
CopyToHost(Averages.Phase_tplus.get(),Phi,N*sizeof(double));
//...........................................................................
//...........................................................................
// Copy the data for for the analysis timestep
@ -1136,11 +1136,11 @@ int main(int argc, char **argv)
//...........................................................................
DeviceBarrier();
ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
CopyToHost(Averages.Phase.data,Phi,N*sizeof(double));
CopyToHost(Averages.Press.data,Pressure,N*sizeof(double));
CopyToHost(Averages.Vel_x.data,&Velocity[0],N*sizeof(double));
CopyToHost(Averages.Vel_y.data,&Velocity[N],N*sizeof(double));
CopyToHost(Averages.Vel_z.data,&Velocity[2*N],N*sizeof(double));
CopyToHost(Averages.Phase.get(),Phi,N*sizeof(double));
CopyToHost(Averages.Press.get(),Pressure,N*sizeof(double));
CopyToHost(Averages.Vel_x.get(),&Velocity[0],N*sizeof(double));
CopyToHost(Averages.Vel_y.get(),&Velocity[N],N*sizeof(double));
CopyToHost(Averages.Vel_z.get(),&Velocity[2*N],N*sizeof(double));
MPI_Barrier(MPI_COMM_WORLD);
//...........................................................................
@ -1551,8 +1551,8 @@ int main(int argc, char **argv)
//...........................................................................
// Copy the phase indicator field for the earlier timestep
DeviceBarrier();
CopyToHost(Averages.Phase.data,Phi,N*sizeof(double));
Averages.ColorToSignedDistance(beta,Averages.Phase.data,Averages.Phase_tplus.data);
CopyToHost(Averages.Phase.get(),Phi,N*sizeof(double));
Averages.ColorToSignedDistance(beta,Averages.Phase.get(),Averages.Phase_tplus.get());
//...........................................................................
}
if (timestep%5000 == 0){
@ -1563,23 +1563,23 @@ int main(int argc, char **argv)
//...........................................................................
DeviceBarrier();
ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
CopyToHost(Averages.Phase.data,Phi,N*sizeof(double));
CopyToHost(Averages.Press.data,Pressure,N*sizeof(double));
CopyToHost(Averages.Vel_x.data,&Velocity[0],N*sizeof(double));
CopyToHost(Averages.Vel_y.data,&Velocity[N],N*sizeof(double));
CopyToHost(Averages.Vel_z.data,&Velocity[2*N],N*sizeof(double));
CopyToHost(Averages.Phase.get(),Phi,N*sizeof(double));
CopyToHost(Averages.Press.get(),Pressure,N*sizeof(double));
CopyToHost(Averages.Vel_x.get(),&Velocity[0],N*sizeof(double));
CopyToHost(Averages.Vel_y.get(),&Velocity[N],N*sizeof(double));
CopyToHost(Averages.Vel_z.get(),&Velocity[2*N],N*sizeof(double));
MPI_Barrier(MPI_COMM_WORLD);
}
if (timestep%5000 == 5){
//...........................................................................
// Copy the phase indicator field for the later timestep
DeviceBarrier();
CopyToHost(Averages.Phase_tminus.data,Phi,N*sizeof(double));
Averages.ColorToSignedDistance(beta,Averages.Phase_tminus.data,Averages.Phase_tminus.data);
CopyToHost(Averages.Phase_tminus.get(),Phi,N*sizeof(double));
Averages.ColorToSignedDistance(beta,Averages.Phase_tminus.get(),Averages.Phase_tminus.get());
//....................................................................
Averages.Initialize();
Averages.ComputeDelPhi();
Averages.ColorToSignedDistance(beta,Averages.Phase.data,Averages.SDn.data);
Averages.ColorToSignedDistance(beta,Averages.Phase.get(),Averages.SDn.get());
Averages.UpdateMeshValues();
Averages.ComputeLocal();
Averages.Reduce();
@ -1747,22 +1747,22 @@ int main(int argc, char **argv)
sprintf(LocalRankFilename,"%s%s","Phase.",LocalRankString);
FILE *PHASE;
PHASE = fopen(LocalRankFilename,"wb");
fwrite(Averages.Phase.data,8,N,PHASE);
// fwrite(MeanCurvature.data,8,N,PHASE);
fwrite(Averages.Phase.get(),8,N,PHASE);
// fwrite(MeanCurvature.get(),8,N,PHASE);
fclose(PHASE);
//#endif
sprintf(LocalRankFilename,"%s%s","Pressure.",LocalRankString);
FILE *PRESS;
PRESS = fopen(LocalRankFilename,"wb");
fwrite(Averages.Press.data,8,N,PRESS);
fwrite(Averages.Press.get(),8,N,PRESS);
fclose(PRESS);
/* sprintf(LocalRankFilename,"%s%s","dPdt.",LocalRankString);
FILE *SPEED;
SPEED = fopen(LocalRankFilename,"wb");
fwrite(dPdt.data,8,N,SPEED);
fwrite(dPdt.get(),8,N,SPEED);
fclose(SPEED);
sprintf(LocalRankFilename,"%s%s","DenA.",LocalRankString);
@ -1787,7 +1787,7 @@ int main(int argc, char **argv)
}
}
}
fwrite(Phase.data,8,N,GRAD);
fwrite(Phase.get(),8,N,GRAD);
fclose(GRAD);
*/
/* double *DensityValues;

8
tests/lbpm_disc_pp.cpp
View File

@ -287,7 +287,7 @@ int main(int argc, char **argv)
}
//.......................................................................
SignedDistanceDiscPack(SignDist.data,ndiscs,cx,cy,rad,Lx,Ly,Lz,Nx,Ny,Nz,
SignedDistanceDiscPack(SignDist.get(),ndiscs,cx,cy,rad,Lx,Ly,Lz,Nx,Ny,Nz,
iproc,jproc,kproc,nprocx,nprocy,nprocz);
//.......................................................................
// Assign walls in the signed distance functions (x,y boundaries)
@ -325,11 +325,11 @@ int main(int argc, char **argv)
for ( j=1;j<Ny-1;j++){
for ( i=1;i<Nx-1;i++){
n = k*Nx*Ny+j*Nx+i;
if (SignDist.data[n] > 0.0){
if (SignDist(n) > 0.0){
id[n] = 2;
}
// compute the porosity (actual interface location used)
if (SignDist.data[n] > 0.0){
if (SignDist(n) > 0.0){
sum++;
}
}
@ -363,7 +363,7 @@ int main(int argc, char **argv)
//.......................................................................
sprintf(LocalRankString,"%05d",rank);
sprintf(LocalRankFilename,"%s%s","SignDist.",LocalRankString);
WriteLocalSolidDistance(LocalRankFilename, SignDist.data, N);
WriteLocalSolidDistance(LocalRankFilename, SignDist.get(), N);
//......................................................................
// ****************************************************

8
tests/lbpm_sphere_pp.cpp
View File

@ -185,7 +185,7 @@ int main(int argc, char **argv)
MPI_Bcast(&D,1,MPI_DOUBLE,0,MPI_COMM_WORLD);
//.......................................................................
SignedDistance(SignDist.data,nspheres,cx,cy,cz,rad,Lx,Ly,Lz,Nx,Ny,Nz,
SignedDistance(SignDist.get(),nspheres,cx,cy,cz,rad,Lx,Ly,Lz,Nx,Ny,Nz,
iproc,jproc,kproc,nprocx,nprocy,nprocz);
//.......................................................................
// Assign the phase ID field based on the signed distance
@ -204,11 +204,11 @@ int main(int argc, char **argv)
for ( j=1;j<Ny-1;j++){
for ( i=1;i<Nx-1;i++){
n = k*Nx*Ny+j*Nx+i;
if (SignDist.data[n] > 0.0){
if (SignDist(n) > 0.0){
id[n] = 2;
}
// compute the porosity (actual interface location used)
if (SignDist.data[n] > 0.0){
if (SignDist(n) > 0.0){
sum++;
}
}
@ -242,7 +242,7 @@ int main(int argc, char **argv)
//.......................................................................
sprintf(LocalRankString,"%05d",rank);
sprintf(LocalRankFilename,"%s%s","SignDist.",LocalRankString);
WriteLocalSolidDistance(LocalRankFilename, SignDist.data, N);
WriteLocalSolidDistance(LocalRankFilename, SignDist.get(), N);
//......................................................................
// ****************************************************