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

Moving to a new signed distance function

Browse Source
This commit is contained in:
Mark Berrill
2018-05-17 12:07:20 -04:00
parent 5a686278bf
commit 195481f615
20 changed files with 502 additions and 659 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
analysis/TwoPhase.cpp
View File

@@ -232,7 +232,7 @@ void TwoPhase::ColorToSignedDistance(double Beta, DoubleArray &ColorData, Double
}
}
Eikonal(DistData,TempID,Dm,50);
CalcDist(DistData,TempID,Dm);
for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){

4
analysis/analysis.cpp
View File

@@ -329,7 +329,7 @@ static int LocalToGlobalIDs( int nx, int ny, int nz, const RankInfoStruct& rank_
}
const BlobIDArray LocalIDs = IDs;
// Copy the ids and get the neighbors through the halos
fillHalo<BlobIDType> fillData(comm,rank_info,nx,ny,nz,1,1,1,0,1,true,true,true);
fillHalo<BlobIDType> fillData(comm,rank_info,{nx,ny,nz},{1,1,1},0,1,{true,true,true});
fillData.fill(IDs);
// Create a list of all neighbor ranks (excluding self)
std::vector<int> neighbors;
@@ -420,7 +420,7 @@ static int LocalToGlobalIDs( int nx, int ny, int nz, const RankInfoStruct& rank_
}
}
// Fill the ghosts
fillHalo<BlobIDType> fillData2(comm,rank_info,nx,ny,nz,1,1,1,0,1,true,true,true);
fillHalo<BlobIDType> fillData2(comm,rank_info,{nx,ny,nz},{1,1,1},0,1,{true,true,true});
fillData2.fill(IDs);
// Reorder based on size (and compress the id space
int N_blobs_global = ReorderBlobIDs2(IDs,N_blobs_tot,ngx,ngy,ngz,comm);

213
analysis/distance.cpp Normal file
View File

@@ -0,0 +1,213 @@
#include "analysis/distance.h"
// Check if we need to recompute distance after updating ghose values
static bool checkUpdate( const Array<Vec> &d, double dx, double dy, double dz )
{
auto s = d.size();
bool test[3] = { false, false, false };
// Check x-direction
Vec v1, v2;
for (size_t k=1; k<s[2]-1; k++) {
for (size_t j=1; j<s[1]-1; j++) {
v1 = d(1,j,k);
v2 = d(0,j,k);
v2.x += dx;
test[0] = test[0] || v2 < v1;
v1 = d(s[0]-2,j,k);
v2 = d(s[0]-1,j,k);
v2.x -= dx;
test[0] = test[0] || v2 < v1;
}
}
// Check y-direction
for (size_t k=1; k<s[2]-1; k++) {
for (size_t i=1; i<s[0]-1; i++) {
v1 = d(i,1,k);
v2 = d(i,0,k);
v2.y += dy;
test[1] = test[1] || v2 < v1;
v1 = d(i,s[1]-2,k);
v2 = d(i,s[1]-1,k);
v2.y -= dy;
test[1] = test[1] || v2 < v1;
}
}
// Check z-direction
for (size_t j=1; j<s[1]-1; j++) {
for (size_t i=1; i<s[0]-1; i++) {
v1 = d(i,j,1);
v2 = d(i,j,0);
v2.z += dz;
test[2] = test[2] || v2 < v1;
v1 = d(i,j,s[2]-2);
v2 = d(i,j,s[2]-1);
v2.z -= dz;
test[2] = test[2] || v2 < v1;
}
}
return test[0] || test[1] || test[2];
}
/******************************************************************
* A fast distance calculation *
******************************************************************/
template<class TYPE>
void CalcDist( Array<TYPE> &Distance, const Array<char> &ID, const Domain &Dm, const std::array<bool,3>& periodic )
{
ASSERT( Distance.size() == ID.size() );
std::array<int,3> n = { Dm.Nx-2, Dm.Ny-2, Dm.Nz-2 };
fillHalo<int> fillData( Dm.Comm, Dm.rank_info, n, {1,1,1}, 50, 1, {true,false,false}, periodic );
Array<int> id(ID.size());
Array<Vec> vecDist(Distance.size());
for (size_t i=0; i<ID.length(); i++)
id(i) = ID(i) == 0 ? -1:1;
fillData.fill( id );
CalcVecDist( vecDist, id, Dm, periodic );
for (size_t i=0; i<Distance.length(); i++)
Distance(i) = id(i)*vecDist(i).norm();
}
void CalcVecDist( Array<Vec> &d, const Array<int> &ID0, const Domain &Dm, const std::array<bool,3>& periodic )
{
const double dx = 1.0;
const double dy = 1.0;
const double dz = 1.0;
std::array<int,3> N = { Dm.Nx, Dm.Ny, Dm.Nz };
std::array<int,3> n = { Dm.Nx-2, Dm.Ny-2, Dm.Nz-2 };
// Create ID with ghosts
Array<int> ID(N[0],N[1],N[2]);
fillHalo<int> fillDataID( Dm.Comm, Dm.rank_info, n, {1,1,1}, 50, 1, {true,true,true}, periodic );
fillDataID.copy( ID0, ID );
// Fill ghosts with nearest neighbor
for (int k=1; k<N[2]-1; k++) {
for (int j=1; j<N[1]-1; j++) {
ID(0,j,k) = ID(1,j,k);
ID(N[0]-1,j,k) = ID(N[0]-2,j,k);
}
}
for (int k=1; k<N[2]-1; k++) {
for (int i=0; i<N[0]; i++) {
ID(i,0,k) = ID(i,1,k);
ID(i,N[1]-1,k) = ID(i,N[1]-2,k);
}
}
for (int i=0; i<N[0]; i++) {
for (int j=0; j<N[1]; j++) {
ID(i,j,0) = ID(i,j,1);
ID(i,j,N[1]-1) = ID(i,j,N[2]-2);
}
}
// Communicate ghosts
fillDataID.fill( ID );
// Create communicator for distance
fillHalo<Vec> fillData( Dm.Comm, Dm.rank_info, n, {1,1,1}, 50, 1, {true,false,false}, periodic );
// Initialize the vector distance
d.fill( Vec( 1e50, 1e50, 1e50 ) );
const double dx0 = 0.5*dx;
const double dy0 = 0.5*dy;
const double dz0 = 0.5*dz;
//const double dxy0 = 0.25*sqrt( dx*dx + dy*dy );
//const double dxz0 = 0.25*sqrt( dx*dx + dz*dz );
//const double dyz0 = 0.25*sqrt( dy*dy + dz*dz );
//const double dxyz0 = sqrt( dx*dx + dy*dy + dz*dz );
for (int k=1; k<Dm.Nz-1; k++) {
for (int j=1; j<Dm.Ny-1; j++) {
for (int i=1; i<Dm.Nx-1; i++) {
int id = ID(i,j,k);
bool x[2] = { id != ID(i-1,j,k), id != ID(i+1,j,k) };
bool y[2] = { id != ID(i,j-1,k), id != ID(i,j+1,k) };
bool z[2] = { id != ID(i,j,k-1), id != ID(i,j,k+1) };
if ( x[0] ) d(i,j,k) = Vec( dx0, 0, 0 );
if ( x[1] ) d(i,j,k) = Vec( -dx0, 0, 0 );
if ( y[0] ) d(i,j,k) = Vec( 0, dy0, 0 );
if ( y[1] ) d(i,j,k) = Vec( 0, -dy0, 0 );
if ( z[0] ) d(i,j,k) = Vec( 0, 0, dz0 );
if ( z[1] ) d(i,j,k) = Vec( 0, 0, -dz0 );
/*if ( x[0] && y[0] ) d(i,j,k) = Vec( dxy0, dxy0, 0 );
if ( x[0] && y[1] ) d(i,j,k) = Vec( dxy0, -dxy0, 0 );
if ( x[1] && y[0] ) d(i,j,k) = Vec( -dxy0, dxy0, 0 );
if ( x[1] && y[1] ) d(i,j,k) = Vec( -dxy0, -dxy0, 0 );
if ( x[0] && z[0] ) d(i,j,k) = Vec( dxz0, 0, dxz0 );
if ( x[0] && z[1] ) d(i,j,k) = Vec( dxz0, 0, -dxz0 );
if ( x[1] && z[0] ) d(i,j,k) = Vec( -dxz0, 0, dxz0 );
if ( x[1] && z[1] ) d(i,j,k) = Vec( -dxz0, 0, -dxz0 );
if ( y[0] && z[0] ) d(i,j,k) = Vec( 0, dyz0, dyz0 );
if ( y[0] && z[1] ) d(i,j,k) = Vec( 0, dyz0, -dyz0 );
if ( y[1] && z[0] ) d(i,j,k) = Vec( 0, -dyz0, dyz0 );
if ( y[1] && z[1] ) d(i,j,k) = Vec( 0, -dyz0, -dyz0 );*/
}
}
}
int N_it = Dm.nprocx() + Dm.nprocy() + Dm.nprocz() + 3;
for ( int it=0; it<N_it; it++ ) {
//if ( Dm.rank() == 0 )
// printf("Runing iteration %i\n",it+1);
// Propagate +/- x-direction
for (int k=0; k<Dm.Nz; k++) {
for (int j=0; j<Dm.Ny; j++) {
for (int i=1; i<Dm.Nx; i++) {
auto v1 = d(i,j,k);
auto v2 = d(i-1,j,k);
v2.x += dx;
if ( v2 < v1 )
d(i,j,k) = v2;
}
for (int i=Dm.Nx-2; i>=0; i--) {
auto v1 = d(i,j,k);
auto v2 = d(i+1,j,k);
v2.x -= dx;
if ( v2 < v1 )
d(i,j,k) = v2;
}
}
}
// Propagate +/- y-direction
for (int k=0; k<Dm.Nz; k++) {
for (int i=0; i<Dm.Nx; i++) {
for (int j=1; j<Dm.Ny; j++) {
auto v1 = d(i,j,k);
auto v2 = d(i,j-1,k);
v2.y += dy;
if ( v2 < v1 )
d(i,j,k) = v2;
}
for (int j=Dm.Ny-2; j>=0; j--) {
auto v1 = d(i,j,k);
auto v2 = d(i,j+1,k);
v2.y -= dy;
if ( v2 < v1 )
d(i,j,k) = v2;
}
}
}
// Propagate +/- z-direction
for (int j=0; j<Dm.Ny; j++) {
for (int i=0; i<Dm.Nx; i++) {
for (int k=1; k<Dm.Nz; k++) {
auto v1 = d(i,j,k);
auto v2 = d(i,j,k-1);
v2.z += dz;
if ( v2 < v1 )
d(i,j,k) = v2;
}
for (int k=Dm.Nz-2; k>=0; k--) {
auto v1 = d(i,j,k);
auto v2 = d(i,j,k+1);
v2.z -= dz;
if ( v2 < v1 )
d(i,j,k) = v2;
}
}
}
fillData.fill( d );
bool test = checkUpdate( d, dx, dy, dz );
test = sumReduce( Dm.Comm, test );
if ( !test )
break;
}
}
template void CalcDist<float>( Array<float>&, const Array<char>&, const Domain&, const std::array<bool,3>& );
template void CalcDist<double>( Array<double>&, const Array<char>&, const Domain&, const std::array<bool,3>& );

40
analysis/distance.h Normal file
View File

@@ -0,0 +1,40 @@
#ifndef Distance_H_INC
#define Distance_H_INC
#include "common/Domain.h"
struct Vec {
double x;
double y;
double z;
inline Vec(): x(0), y(0), z(0) {}
inline Vec( double x_, double y_, double z_ ): x(x_), y(y_), z(z_) {}
inline double norm() const { return sqrt(x*x+y*y+z*z); }
};
inline bool operator<(const Vec& l, const Vec& r){ return l.x*l.x+l.y*l.y+l.z*l.z < r.x*r.x+r.y*r.y+r.z*r.z; }
/*!
* @brief Calculate the distance using a simple method
* @details This routine calculates the vector distance to the nearest domain surface.
* @param[out] Distance Distance function
* @param[in] ID Segmentation id
* @param[in] Dm Domain information
* @param[in] periodic Directions that are periodic
*/
template<class TYPE>
void CalcDist( Array<TYPE> &Distance, const Array<char> &ID, const Domain &Dm, const std::array<bool,3>& periodic = {true,true,true} );
/*!
* @brief Calculate the distance using a simple method
* @details This routine calculates the vector distance to the nearest domain surface.
* @param[out] Distance Distance function
* @param[in] ID Domain id
* @param[in] Dm Domain information
* @param[in] periodic Directions that are periodic
*/
void CalcVecDist( Array<Vec> &Distance, const Array<int> &ID, const Domain &Dm, const std::array<bool,3>& periodic );
#endif

384
analysis/eikonal.cpp
View File

@@ -1,384 +0,0 @@
#include "analysis/eikonal.h"
#include "analysis/imfilter.h"
static inline float minmod(float &a, float &b)
{
float value = a;
if ( a*b < 0.0)
value=0.0;
else if (fabs(a) > fabs(b))
value = b;
return value;
}
static inline double minmod(double &a, double &b){
double value;
value = a;
if ( a*b < 0.0) value=0.0;
else if (fabs(a) > fabs(b)) value = b;
return value;
}
template<class TYPE>
static inline MPI_Datatype getType( );
template<> inline MPI_Datatype getType<float>() { return MPI_FLOAT; }
template<> inline MPI_Datatype getType<double>() { return MPI_DOUBLE; }
/******************************************************************
* Solve the eikonal equation *
******************************************************************/
template<class TYPE>
TYPE Eikonal( Array<TYPE> &Distance, const Array<char> &ID, const Domain &Dm, const int timesteps)
{
/*
* This routine converts the data in the Distance array to a signed distance
* by solving the equation df/dt = sign(1-|grad f|), where Distance provides
* the values of f on the mesh associated with domain Dm
* It has been tested with segmented data initialized to values [-1,1]
* and will converge toward the signed distance to the surface bounding the associated phases
*
* Reference:
* Min C (2010) On reinitializing level set functions, Journal of Computational Physics 229
*/
int xdim = Dm.Nx-2;
int ydim = Dm.Ny-2;
int zdim = Dm.Nz-2;
fillHalo<TYPE> fillData(Dm.Comm, Dm.rank_info,xdim,ydim,zdim,1,1,1,0,1);
// Arrays to store the second derivatives
Array<TYPE> Dxx(Dm.Nx,Dm.Ny,Dm.Nz);
Array<TYPE> Dyy(Dm.Nx,Dm.Ny,Dm.Nz);
Array<TYPE> Dzz(Dm.Nx,Dm.Ny,Dm.Nz);
Array<int8_t> sign(ID.size());
for (size_t i=0; i<sign.length(); i++)
sign(i) = ID(i) == 0 ? -1:1;
int count = 0;
double dt = 0.1;
double GlobalVar = 0.0;
while (count < timesteps){
// Communicate the halo of values
fillData.fill(Distance);
// Compute second order derivatives
for (int k=1;k<Dm.Nz-1;k++){
for (int j=1;j<Dm.Ny-1;j++){
for (int i=1;i<Dm.Nx-1;i++){
Dxx(i,j,k) = Distance(i+1,j,k) + Distance(i-1,j,k) - 2*Distance(i,j,k);
Dyy(i,j,k) = Distance(i,j+1,k) + Distance(i,j-1,k) - 2*Distance(i,j,k);
Dzz(i,j,k) = Distance(i,j,k+1) + Distance(i,j,k-1) - 2*Distance(i,j,k);
}
}
}
fillData.fill(Dxx);
fillData.fill(Dyy);
fillData.fill(Dzz);
double LocalMax = 0.0;
double LocalVar = 0.0;
// Execute the next timestep
for (int k=1;k<Dm.Nz-1;k++){
for (int j=1;j<Dm.Ny-1;j++){
for (int i=1;i<Dm.Nx-1;i++){
double s = sign(i,j,k);
// local second derivative terms
double Dxxp = minmod(Dxx(i,j,k),Dxx(i+1,j,k));
double Dyyp = minmod(Dyy(i,j,k),Dyy(i,j+1,k));
double Dzzp = minmod(Dzz(i,j,k),Dzz(i,j,k+1));
double Dxxm = minmod(Dxx(i,j,k),Dxx(i-1,j,k));
double Dyym = minmod(Dyy(i,j,k),Dyy(i,j-1,k));
double Dzzm = minmod(Dzz(i,j,k),Dzz(i,j,k-1));
/* //............Compute upwind derivatives ...................
Dxp = Distance(i+1,j,k) - Distance(i,j,k) + 0.5*Dxxp;
Dyp = Distance(i,j+1,k) - Distance(i,j,k) + 0.5*Dyyp;
Dzp = Distance(i,j,k+1) - Distance(i,j,k) + 0.5*Dzzp;
Dxm = Distance(i,j,k) - Distance(i-1,j,k) + 0.5*Dxxm;
Dym = Distance(i,j,k) - Distance(i,j-1,k) + 0.5*Dyym;
Dzm = Distance(i,j,k) - Distance(i,j,k-1) + 0.5*Dzzm;
*/
double Dxp = Distance(i+1,j,k)- Distance(i,j,k) - 0.5*Dxxp;
double Dyp = Distance(i,j+1,k)- Distance(i,j,k) - 0.5*Dyyp;
double Dzp = Distance(i,j,k+1)- Distance(i,j,k) - 0.5*Dzzp;
double Dxm = Distance(i,j,k) - Distance(i-1,j,k) + 0.5*Dxxm;
double Dym = Distance(i,j,k) - Distance(i,j-1,k) + 0.5*Dyym;
double Dzm = Distance(i,j,k) - Distance(i,j,k-1) + 0.5*Dzzm;
// Compute upwind derivatives for Godunov Hamiltonian
double Dx, Dy, Dz;
if ( s < 0.0){
if (Dxp + Dxm > 0.f)
Dx = Dxp*Dxp;
else
Dx = Dxm*Dxm;
if (Dyp + Dym > 0.f)
Dy = Dyp*Dyp;
else
Dy = Dym*Dym;
if (Dzp + Dzm > 0.f)
Dz = Dzp*Dzp;
else
Dz = Dzm*Dzm;
}
else{
if (Dxp + Dxm < 0.f)
Dx = Dxp*Dxp;
else
Dx = Dxm*Dxm;
if (Dyp + Dym < 0.f)
Dy = Dyp*Dyp;
else
Dy = Dym*Dym;
if (Dzp + Dzm < 0.f)
Dz = Dzp*Dzp;
else
Dz = Dzm*Dzm;
}
//Dx = max(Dxp*Dxp,Dxm*Dxm);
//Dy = max(Dyp*Dyp,Dym*Dym);
//Dz = max(Dzp*Dzp,Dzm*Dzm);
double norm = sqrt(Dx + Dy + Dz);
if (norm > 1.0)
norm = 1.0;
Distance(i,j,k) += dt*s*(1.0 - norm);
LocalVar += dt*s*(1.0 - norm);
if (fabs(dt*s*(1.0 - norm)) > LocalMax)
LocalMax = fabs(dt*s*(1.0 - norm));
}
}
}
double GlobalMax;
MPI_Allreduce(&LocalVar,&GlobalVar,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&LocalMax,&GlobalMax,1,MPI_DOUBLE,MPI_MAX,Dm.Comm);
GlobalVar /= (Dm.Nx-2)*(Dm.Ny-2)*(Dm.Nz-2)*Dm.nprocx()*Dm.nprocy()*Dm.nprocz();
count++;
if (count%50 == 0 && Dm.rank()==0 ){
printf("Time=%i, Max variation=%f, Global variation=%f \n",count,GlobalMax,GlobalVar);
fflush(stdout);
}
if (fabs(GlobalMax) < 1e-5){
if (Dm.rank()==0)
printf("Exiting with max tolerance of 1e-5 \n");
count=timesteps;
}
}
return GlobalVar;
}
template float Eikonal<float>( Array<float>&, const Array<char>&, const Domain&, int );
template double Eikonal<double>( Array<double>&, const Array<char>&, const Domain&, int );
/******************************************************************
* A fast distance calculation *
******************************************************************/
bool CalcDist3DIteration( Array<float> &Distance, const Domain & )
{
const float sq2 = sqrt(2.0f);
const float sq3 = sqrt(3.0f);
float dist0[27];
dist0[0] = sq3; dist0[1] = sq2; dist0[2] = sq3;
dist0[3] = sq2; dist0[4] = 1; dist0[5] = sq2;
dist0[6] = sq3; dist0[7] = sq2; dist0[8] = sq3;
dist0[9] = sq2; dist0[10] = 1; dist0[11] = sq2;
dist0[12] = 1; dist0[13] = 0; dist0[14] = 1;
dist0[15] = sq2; dist0[16] = 1; dist0[17] = sq2;
dist0[18] = sq3; dist0[19] = sq2; dist0[20] = sq3;
dist0[21] = sq2; dist0[22] = 1; dist0[23] = sq2;
dist0[24] = sq3; dist0[25] = sq2; dist0[26] = sq3;
bool changed = false;
for (size_t k=1; k<Distance.size(2)-1; k++) {
for (size_t j=1; j<Distance.size(1)-1; j++) {
for (size_t i=1; i<Distance.size(0)-1; i++) {
float dist[27];
dist[0] = Distance(i-1,j-1,k-1); dist[1] = Distance(i,j-1,k-1); dist[2] = Distance(i+1,j-1,k-1);
dist[3] = Distance(i-1,j,k-1); dist[4] = Distance(i,j,k-1); dist[5] = Distance(i+1,j,k-1);
dist[6] = Distance(i-1,j+1,k-1); dist[7] = Distance(i,j+1,k-1); dist[8] = Distance(i+1,j+1,k-1);
dist[9] = Distance(i-1,j-1,k); dist[10] = Distance(i,j-1,k); dist[11] = Distance(i+1,j-1,k);
dist[12] = Distance(i-1,j,k); dist[13] = Distance(i,j,k); dist[14] = Distance(i+1,j,k);
dist[15] = Distance(i-1,j+1,k); dist[16] = Distance(i,j+1,k); dist[17] = Distance(i+1,j+1,k);
dist[18] = Distance(i-1,j-1,k+1); dist[19] = Distance(i,j-1,k+1); dist[20] = Distance(i+1,j-1,k+1);
dist[21] = Distance(i-1,j,k+1); dist[22] = Distance(i,j,k+1); dist[23] = Distance(i+1,j,k+1);
dist[24] = Distance(i-1,j+1,k+1); dist[25] = Distance(i,j+1,k+1); dist[26] = Distance(i+1,j+1,k+1);
float tmp = 1e100;
for (int k=0; k<27; k++)
tmp = std::min(tmp,dist[k]+dist0[k]);
if ( tmp < Distance(i,j,k) ) {
Distance(i,j,k) = tmp;
changed = true;
}
}
}
}
return changed;
}
void CalcDist3D( Array<float> &Distance, const Array<char> &ID, const Domain &Dm )
{
PROFILE_START("Calc Distance");
// Initialize the distance to be 0 fore the cells adjacent to the interface
Distance.fill(1e100);
for (size_t k=1; k<ID.size(2)-1; k++) {
for (size_t j=1; j<ID.size(1)-1; j++) {
for (size_t i=1; i<ID.size(0)-1; i++) {
char id = ID(i,j,k);
if ( id!=ID(i-1,j,k) || id!=ID(i+1,j,k) || id!=ID(i,j-1,k) || id!=ID(i,j+1,k) || id!=ID(i,j,k-1) || id!=ID(i,j,k+1) )
Distance(i,j,k) = 0.5;
}
}
}
// Compute the distance everywhere
fillHalo<float> fillData(Dm.Comm, Dm.rank_info,Dm.Nx,Dm.Ny,Dm.Nz,1,1,1,0,1);
while ( true ) {
// Communicate the halo of values
fillData.fill(Distance);
// The distance of the cell is the minimum of the distance of the neighbors plus the distance to that node
bool changed = CalcDist3DIteration( Distance, Dm );
changed = sumReduce(Dm.Comm,changed);
if ( !changed )
break;
}
// Update the sign of the distance
for (size_t i=0; i<ID.length(); i++)
Distance(i) *= ID(i)>0 ? 1:-1;
PROFILE_STOP("Calc Distance");
}
/******************************************************************
* A fast distance calculation *
******************************************************************/
void CalcDistMultiLevelHelper( Array<float> &Distance, const Domain &Dm )
{
size_t ratio = 4;
std::function<float(const Array<float>&)> coarsen = [ratio]( const Array<float>& data )
{
float tmp = 1e100;
int nx = data.size(0);
int ny = data.size(1);
int nz = data.size(2);
for (int k=0; k<nz; k++) {
float z = k-0.5*(nz-1);
for (int j=0; j<ny; j++) {
float y = j-0.5*(ny-1);
for (int i=0; i<nx; i++) {
float x = i-0.5*(nx-1);
tmp = std::min(data(i,j,k)+sqrt(x*x+y*y+z*z),tmp);
}
}
}
return tmp/ratio;
};
int Nx = Dm.Nx-2;
int Ny = Dm.Ny-2;
int Nz = Dm.Nz-2;
ASSERT(int(Distance.size(0))==Nx+2&&int(Distance.size(1))==Ny+2&&int(Distance.size(2))==Nz+2);
fillHalo<float> fillData(Dm.Comm,Dm.rank_info,Nx,Ny,Nz,1,1,1,0,1);
if ( Nx%ratio==0 && Nx>8 && Ny%ratio==0 && Ny>8 && Nz%ratio==0 && Nz>8 ) {
// Use recursive version
int Nr = std::max(std::max(ratio,ratio),ratio);
// Run Nr iterations, communicate, run Nr iterations
for (int i=0; i<Nr; i++)
CalcDist3DIteration( Distance, Dm );
/*fillData.fill(Distance);
for (int i=0; i<Nr; i++)
CalcDist3DIteration( Distance, Dm );*/
// Coarsen
Array<float> dist(Nx,Ny,Nz);
fillData.copy(Distance,dist);
auto db = Dm.getDatabase()->cloneDatabase();
auto n = db->getVector<int>( "n" );
db->putVector<int>( "n", { n[0]/ratio, n[1]/ratio, n[2]/ratio } );
Domain Dm2(db);
Dm2.CommInit(Dm.Comm);
fillHalo<float> fillData2(Dm2.Comm,Dm2.rank_info,Nx/ratio,Ny/ratio,Nz/ratio,1,1,1,0,1);
auto dist2 = dist.coarsen( {ratio,ratio,ratio}, coarsen );
Array<float> Distance2(Nx/ratio+2,Ny/ratio+2,Nz/ratio+2);
fillData2.copy(dist2,Distance2);
// Solve
CalcDistMultiLevelHelper( Distance2, Dm2 );
// Interpolate the coarse grid to the fine grid
fillData2.copy(Distance2,dist2);
for (int k=0; k<Nz; k++) {
int k2 = k/ratio;
float z = (k-k2*ratio)-0.5*(ratio-1);
for (int j=0; j<Ny; j++) {
int j2 = j/ratio;
float y = (j-j2*ratio)-0.5*(ratio-1);
for (int i=0; i<Nx; i++) {
int i2 = i/ratio;
float x = (i-i2*ratio)-0.5*(ratio-1);
dist(i,j,k) = std::min(dist(i,j,k),ratio*dist2(i2,j2,k2)+sqrt(x*x+y*y+z*z));
}
}
}
fillData.copy(dist,Distance);
// Run Nr iterations, communicate, run Nr iterations
for (int i=0; i<Nr; i++)
CalcDist3DIteration( Distance, Dm );
fillData.fill(Distance);
for (int i=0; i<Nr; i++)
CalcDist3DIteration( Distance, Dm );
} else {
// Use coarse-grid version
while ( true ) {
// Communicate the halo of values
fillData.fill(Distance);
// The distance of the cell is the minimum of the distance of the neighbors plus the distance to that node
bool changed = CalcDist3DIteration( Distance, Dm );
changed = sumReduce(Dm.Comm,changed);
if ( !changed )
break;
}
}
}
void CalcDistMultiLevel( Array<float> &Distance, const Array<char> &ID, const Domain &Dm )
{
PROFILE_START("Calc Distance Multilevel");
int Nx = Dm.Nx-2;
int Ny = Dm.Ny-2;
int Nz = Dm.Nz-2;
ASSERT(int(Distance.size(0))==Nx+2&&int(Distance.size(1))==Ny+2&&int(Distance.size(2))==Nz+2);
fillHalo<float> fillData(Dm.Comm,Dm.rank_info,Nx,Ny,Nz,1,1,1,0,1);
// Initialize the distance to be 0 fore the cells adjacent to the interface
Distance.fill(1e100);
for (size_t k=1; k<ID.size(2)-1; k++) {
for (size_t j=1; j<ID.size(1)-1; j++) {
for (size_t i=1; i<ID.size(0)-1; i++) {
char id = ID(i,j,k);
if ( id!=ID(i-1,j,k) || id!=ID(i+1,j,k) || id!=ID(i,j-1,k) || id!=ID(i,j+1,k) || id!=ID(i,j,k-1) || id!=ID(i,j,k+1) )
Distance(i,j,k) = 0.5;
}
}
}
// Solve the for the distance using a recursive method
CalcDistMultiLevelHelper( Distance, Dm );
// Update the sign of the distance
for (size_t i=0; i<ID.length(); i++)
Distance(i) *= ID(i)>0 ? 1:-1;
fillData.fill(Distance);
// Run a quick filter to smooth the data
float sigma = 0.6;
Array<float> H = imfilter::create_filter<float>( { 1 }, "gaussian", &sigma );
std::vector<imfilter::BC> BC(3,imfilter::BC::replicate);
Distance = imfilter::imfilter_separable<float>( Distance, {H,H,H}, BC );
PROFILE_STOP("Calc Distance Multilevel");
}

53
analysis/eikonal.h
View File

@@ -1,53 +0,0 @@
#ifndef Eikonal_H_INC
#define Eikonal_H_INC
#include "common/Domain.h"
/*!
* @brief Calculate the distance solving the Eikonal equation
* @details This routine converts the data in the Distance array to a signed distance
* by solving the equation df/dt = sign*(1-|grad f|), where Distance provides
* the values of f on the mesh associated with domain Dm
* It has been tested with segmented data initialized to values [-1,1]
* and will converge toward the signed distance to the surface bounding the associated phases
*
* Reference:
* Min C (2010) On reinitializing level set functions, Journal of Computational Physics 229
*
* @param[in/out] Distance Distance function
* @param[in] ID Segmentation id
* @param[in] DM Domain information
* @param[in] timesteps Maximum number of timesteps to process
* @return Returns the global variation
*/
template<class TYPE>
TYPE Eikonal( Array<TYPE> &Distance, const Array<char> &ID, const Domain &Dm, int timesteps);
/*!
* @brief Calculate the distance using a simple method
* @details This routine calculates the distance using a very simple method working off the segmentation id.
*
* @param[in/out] Distance Distance function
* @param[in] ID Segmentation id
* @param[in] DM Domain information
* @return Returns the global variation
*/
void CalcDist3D( Array<float> &Distance, const Array<char> &ID, const Domain &Dm );
/*!
* @brief Calculate the distance using a multi-level method
* @details This routine calculates the distance using a multi-grid method
*
* @return Returns the number of cubes in the blob
* @param[in/out] Distance Distance function
* @param[in] ID Segmentation id
* @param[in] DM Domain information
* @return Returns the global variation
*/
void CalcDistMultiLevel( Array<float> &Distance, const Array<char> &ID, const Domain &Dm );
#endif

2
analysis/runAnalysis.cpp
View File

@@ -289,7 +289,7 @@ runAnalysis::runAnalysis( std::shared_ptr<Database> db,
d_ScaLBL_Comm( ScaLBL_Comm ),
d_rank_info( rank_info ),
d_Map( Map ),
d_fillData(Dm.Comm,Dm.rank_info,Dm.Nx-2,Dm.Ny-2,Dm.Nz-2,1,1,1,0,1)
d_fillData(Dm.Comm,Dm.rank_info,{Dm.Nx-2,Dm.Ny-2,Dm.Nz-2},{1,1,1},0,1)
{
NULL_USE( pBC );
INSIST( db, "Input database is empty" );

22
analysis/uCT.cpp
View File

@@ -1,6 +1,6 @@
#include "analysis/uCT.h"
#include "analysis/analysis.h"
#include "analysis/eikonal.h"
#include "analysis/distance.h"
#include "analysis/filters.h"
#include "analysis/imfilter.h"
@@ -149,7 +149,7 @@ void solve( const Array<float>& VOL, Array<float>& Mean, Array<char>& ID,
fillFloat.fill( Mean );
segment( Mean, ID, 0.01 );
// Compute the distance using the segmented volume
Eikonal( Dist, ID, Dm, ID.size(0)*nprocx );
CalcDist( Dist, ID, Dm );
fillFloat.fill(Dist);
smooth( VOL, Dist, 2.0, MultiScaleSmooth, fillFloat );
// Compute non-local mean
@@ -210,9 +210,7 @@ void refine( const Array<float>& Dist_coarse,
//removeDisconnected( ID, Dm );
// Compute the distance using the segmented volume
if ( level > 0 ) {
//Eikonal3D( Dist, ID, Dm, ID.size(0)*nprocx );
//CalcDist3D( Dist, ID, Dm );
CalcDistMultiLevel( Dist, ID, Dm );
CalcDist( Dist, ID, Dm );
fillFloat.fill(Dist);
}
}
@@ -232,7 +230,7 @@ void filter_final( Array<char>& ID, Array<float>& Dist,
int Ny = Dm.Ny-2;
int Nz = Dm.Nz-2;
// Calculate the distance
CalcDistMultiLevel( Dist, ID, Dm );
CalcDist( Dist, ID, Dm );
fillFloat.fill(Dist);
// Compute the range to shrink the volume based on the L2 norm of the distance
Array<float> Dist0(Nx,Ny,Nz);
@@ -256,8 +254,8 @@ void filter_final( Array<char>& ID, Array<float>& Dist,
}
//Array<float> Dist1 = Dist;
//Array<float> Dist2 = Dist;
CalcDistMultiLevel( Dist1, ID1, Dm );
CalcDistMultiLevel( Dist2, ID2, Dm );
CalcDist( Dist1, ID1, Dm );
CalcDist( Dist2, ID2, Dm );
fillFloat.fill(Dist1);
fillFloat.fill(Dist2);
// Keep those regions that are within dx2 of the new volumes
@@ -272,8 +270,8 @@ void filter_final( Array<char>& ID, Array<float>& Dist,
}
}
// Find regions of uncertainty that are entirely contained within another region
fillHalo<double> fillDouble(Dm.Comm,Dm.rank_info,Nx,Ny,Nz,1,1,1,0,1);
fillHalo<BlobIDType> fillInt(Dm.Comm,Dm.rank_info,Nx,Ny,Nz,1,1,1,0,1);
fillHalo<double> fillDouble(Dm.Comm,Dm.rank_info,{Nx,Ny,Nz},{1,1,1},0,1);
fillHalo<BlobIDType> fillInt(Dm.Comm,Dm.rank_info,{Nx,Ny,Nz},{1,1,1},0,1);
BlobIDArray GlobalBlobID;
DoubleArray SignDist(ID.size());
for (size_t i=0; i<ID.length(); i++)
@@ -343,7 +341,7 @@ void filter_final( Array<char>& ID, Array<float>& Dist,
// Perform the final segmentation and update the distance
fillFloat.fill(Mean);
segment( Mean, ID, 0.01 );
CalcDistMultiLevel( Dist, ID, Dm );
CalcDist( Dist, ID, Dm );
fillFloat.fill(Dist);
}
@@ -356,7 +354,7 @@ void filter_src( const Domain& Dm, Array<float>& src )
int Nx = Dm.Nx-2;
int Ny = Dm.Ny-2;
int Nz = Dm.Nz-2;
fillHalo<float> fillFloat(Dm.Comm,Dm.rank_info,Nx,Ny,Nz,1,1,1,0,1);
fillHalo<float> fillFloat(Dm.Comm,Dm.rank_info,{Nx,Ny,Nz},{1,1,1},0,1);
// Perform a hot-spot filter on the data
std::vector<imfilter::BC> BC = { imfilter::BC::replicate, imfilter::BC::replicate, imfilter::BC::replicate };
std::function<float(const Array<float>&)> filter_3D = []( const Array<float>& data )

25
common/Communication.h
View File

@@ -5,6 +5,8 @@
#include "common/Utilities.h"
#include "common/Array.h"
#include <array>
// ********** COMMUNICTION **************************************
/*
//..............................................................
@@ -43,18 +45,17 @@ public:
/*!
* @brief Default constructor
* @param[in] info Rank and neighbor rank info
* @param[in] nx Number of local cells in the x direction
* @param[in] ny Number of local cells in the y direction
* @param[in] nz Number of local cells in the z direction
* @param[in] ngx Number of ghost cells in the x direction
* @param[in] ngy Number of ghost cells in the y direction
* @param[in] ngz Number of ghost cells in the z direction
* @param[in] n Number of local cells
* @param[in] ng Number of ghost cells
* @param[in] tag Initial tag to use for the communication (we will require tag:tag+26)
* @param[in] depth Maximum depth to support
* @param[in] fill Fill {faces,edges,corners}
* @param[in] periodic Periodic dimensions
*/
fillHalo( MPI_Comm comm, const RankInfoStruct& info, int nx, int ny, int nz,
int ngx, int ngy, int ngz, int tag, int depth,
bool fill_face=true, bool fill_edge=true, bool fill_corner=true );
fillHalo( MPI_Comm comm, const RankInfoStruct& info,
std::array<int,3> n, std::array<int,3> ng, int tag, int depth,
std::array<bool,3> fill = {true,true,true},
std::array<bool,3> periodic = {true,true,true} );
//! Destructor
~fillHalo( );
@@ -75,15 +76,17 @@ public:
private:
MPI_Comm comm;
RankInfoStruct info;
int nx, ny, nz, ngx, ngy, ngz, depth;
std::array<int,3> n, ng;
int depth;
bool fill_pattern[3][3][3];
int tag[3][3][3];
int N_send_recv[3][3][3];
TYPE *mem;
TYPE *send[3][3][3], *recv[3][3][3];
MPI_Request send_req[3][3][3], recv_req[3][3][3];
MPI_Comm comm;
size_t N_type;
MPI_Datatype datatype;
fillHalo(); // Private empty constructor
fillHalo(const fillHalo&); // Private copy constructor

157
common/Communication.hpp
View File

@@ -11,16 +11,30 @@
* Structure to store the rank info *
********************************************************/
template<class TYPE>
fillHalo<TYPE>::fillHalo( MPI_Comm comm0, const RankInfoStruct& info0, int nx0, int ny0, int nz0,
int ngx0, int ngy0, int ngz0, int tag0, int depth0,
bool fill_face, bool fill_edge, bool fill_corner ):
info(info0), nx(nx0), ny(ny0), nz(nz0), ngx(ngx0), ngy(ngy0), ngz(ngz0), depth(depth0)
fillHalo<TYPE>::fillHalo( MPI_Comm comm_, const RankInfoStruct& info_,
std::array<int,3> n_, std::array<int,3> ng_, int tag0, int depth_,
std::array<bool,3> fill, std::array<bool,3> periodic ):
comm(comm_), info(info_), n(n_), ng(ng_), depth(depth_)
{
comm = comm0;
datatype = getMPItype<TYPE>();
if ( std::is_same<TYPE,double>() ) {
N_type = 1;
datatype = MPI_DOUBLE;
} else if ( std::is_same<TYPE,float>() ) {
N_type = 1;
datatype = MPI_FLOAT;
} else if ( sizeof(TYPE)%sizeof(double)==0 ) {
N_type = sizeof(TYPE) / sizeof(double);
datatype = MPI_DOUBLE;
} else if ( sizeof(TYPE)%sizeof(float)==0 ) {
N_type = sizeof(TYPE) / sizeof(float);
datatype = MPI_FLOAT;
} else {
N_type = sizeof(TYPE);
datatype = MPI_BYTE;
}
// Set the fill pattern
memset(fill_pattern,0,sizeof(fill_pattern));
if ( fill_face ) {
if ( fill[0] ) {
fill_pattern[0][1][1] = true;
fill_pattern[2][1][1] = true;
fill_pattern[1][0][1] = true;
@@ -28,7 +42,7 @@ fillHalo<TYPE>::fillHalo( MPI_Comm comm0, const RankInfoStruct& info0, int nx0,
fill_pattern[1][1][0] = true;
fill_pattern[1][1][2] = true;
}
if ( fill_edge ) {
if ( fill[1] ) {
fill_pattern[0][0][1] = true;
fill_pattern[0][2][1] = true;
fill_pattern[2][0][1] = true;
@@ -42,7 +56,7 @@ fillHalo<TYPE>::fillHalo( MPI_Comm comm0, const RankInfoStruct& info0, int nx0,
fill_pattern[1][2][0] = true;
fill_pattern[1][2][2] = true;
}
if ( fill_corner ) {
if ( fill[2] ) {
fill_pattern[0][0][0] = true;
fill_pattern[0][0][2] = true;
fill_pattern[0][2][0] = true;
@@ -52,13 +66,50 @@ fillHalo<TYPE>::fillHalo( MPI_Comm comm0, const RankInfoStruct& info0, int nx0,
fill_pattern[2][2][0] = true;
fill_pattern[2][2][2] = true;
}
// Remove communication for non-perioidic directions
if ( !periodic[0] && info.ix==0 ) {
for (int j=0; j<3; j++) {
for (int k=0; k<3; k++)
fill_pattern[0][j][k] = false;
}
}
if ( !periodic[0] && info.ix==info.nx-1 ) {
for (int j=0; j<3; j++) {
for (int k=0; k<3; k++)
fill_pattern[2][j][k] = false;
}
}
if ( !periodic[1] && info.jy==0 ) {
for (int i=0; i<3; i++) {
for (int k=0; k<3; k++)
fill_pattern[i][0][k] = false;
}
}
if ( !periodic[1] && info.jy==info.ny-1 ) {
for (int i=0; i<3; i++) {
for (int k=0; k<3; k++)
fill_pattern[i][2][k] = false;
}
}
if ( !periodic[2] && info.kz==0 ) {
for (int i=0; i<3; i++) {
for (int j=0; j<3; j++)
fill_pattern[i][j][0] = false;
}
}
if ( !periodic[2] && info.kz==info.nz-1 ) {
for (int i=0; i<3; i++) {
for (int j=0; j<3; j++)
fill_pattern[i][j][2] = false;
}
}
// Determine the number of elements for each send/recv
for (int i=0; i<3; i++) {
int ni = (i-1)==0 ? nx:ngx;
int ni = (i-1)==0 ? n[0]:ng[0];
for (int j=0; j<3; j++) {
int nj = (j-1)==0 ? ny:ngy;
int nj = (j-1)==0 ? n[1]:ng[1];
for (int k=0; k<3; k++) {
int nk = (k-1)==0 ? nz:ngz;
int nk = (k-1)==0 ? n[2]:ng[2];
if ( fill_pattern[i][j][k] )
N_send_recv[i][j][k] = ni*nj*nk;
else
@@ -106,9 +157,9 @@ void fillHalo<TYPE>::fill( Array<TYPE>& data )
{
//PROFILE_START("fillHalo::fill",1);
int depth2 = data.size(3);
ASSERT((int)data.size(0)==nx+2*ngx);
ASSERT((int)data.size(1)==ny+2*ngy);
ASSERT((int)data.size(2)==nz+2*ngz);
ASSERT((int)data.size(0)==n[0]+2*ng[0]);
ASSERT((int)data.size(1)==n[1]+2*ng[1]);
ASSERT((int)data.size(2)==n[2]+2*ng[2]);
ASSERT(depth2<=depth);
ASSERT(data.ndim()==3||data.ndim()==4);
// Start the recieves
@@ -117,7 +168,7 @@ void fillHalo<TYPE>::fill( Array<TYPE>& data )
for (int k=0; k<3; k++) {
if ( !fill_pattern[i][j][k] )
continue;
MPI_Irecv( recv[i][j][k], depth2*N_send_recv[i][j][k], datatype,
MPI_Irecv( recv[i][j][k], N_type*depth2*N_send_recv[i][j][k], datatype,
info.rank[i][j][k], tag[2-i][2-j][2-k], comm, &recv_req[i][j][k] );
}
}
@@ -129,7 +180,7 @@ void fillHalo<TYPE>::fill( Array<TYPE>& data )
if ( !fill_pattern[i][j][k] )
continue;
pack( data, i-1, j-1, k-1, send[i][j][k] );
MPI_Isend( send[i][j][k], depth2*N_send_recv[i][j][k], datatype,
MPI_Isend( send[i][j][k], N_type*depth2*N_send_recv[i][j][k], datatype,
info.rank[i][j][k], tag[i][j][k], comm, &send_req[i][j][k] );
}
}
@@ -162,12 +213,12 @@ template<class TYPE>
void fillHalo<TYPE>::pack( const Array<TYPE>& data, int i0, int j0, int k0, TYPE *buffer )
{
int depth2 = data.size(3);
int ni = i0==0 ? nx:ngx;
int nj = j0==0 ? ny:ngy;
int nk = k0==0 ? nz:ngz;
int is = i0==0 ? ngx:((i0==-1)?ngx:nx);
int js = j0==0 ? ngy:((j0==-1)?ngy:ny);
int ks = k0==0 ? ngz:((k0==-1)?ngz:nz);
int ni = i0==0 ? n[0]:ng[0];
int nj = j0==0 ? n[1]:ng[1];
int nk = k0==0 ? n[2]:ng[2];
int is = i0==0 ? ng[0]:((i0==-1)?ng[0]:n[0]);
int js = j0==0 ? ng[1]:((j0==-1)?ng[1]:n[1]);
int ks = k0==0 ? ng[2]:((k0==-1)?ng[2]:n[2]);
for (int d=0; d<depth2; d++) {
for (int k=0; k<nk; k++) {
for (int j=0; j<nj; j++) {
@@ -182,12 +233,12 @@ template<class TYPE>
void fillHalo<TYPE>::unpack( Array<TYPE>& data, int i0, int j0, int k0, const TYPE *buffer )
{
int depth2 = data.size(3);
int ni = i0==0 ? nx:ngx;
int nj = j0==0 ? ny:ngy;
int nk = k0==0 ? nz:ngz;
int is = i0==0 ? ngx:((i0==-1)?0:nx+ngx);
int js = j0==0 ? ngy:((j0==-1)?0:ny+ngy);
int ks = k0==0 ? ngz:((k0==-1)?0:nz+ngz);
int ni = i0==0 ? n[0]:ng[0];
int nj = j0==0 ? n[1]:ng[1];
int nk = k0==0 ? n[2]:ng[2];
int is = i0==0 ? ng[0]:((i0==-1)?0:n[0]+ng[0]);
int js = j0==0 ? ng[1]:((j0==-1)?0:n[1]+ng[1]);
int ks = k0==0 ? ng[2]:((k0==-1)?0:n[2]+ng[2]);
for (int d=0; d<depth2; d++) {
for (int k=0; k<nk; k++) {
for (int j=0; j<nj; j++) {
@@ -208,27 +259,27 @@ template<class TYPE1, class TYPE2>
void fillHalo<TYPE>::copy( const Array<TYPE1>& src, Array<TYPE2>& dst )
{
//PROFILE_START("fillHalo::copy",1);
ASSERT( (int)src.size(0)==nx || (int)src.size(0)==nx+2*ngx );
ASSERT( (int)dst.size(0)==nx || (int)dst.size(0)==nx+2*ngx );
bool src_halo = (int)src.size(0)==nx+2*ngx;
bool dst_halo = (int)dst.size(0)==nx+2*ngx;
ASSERT( (int)src.size(0)==n[0] || (int)src.size(0)==n[0]+2*ng[0] );
ASSERT( (int)dst.size(0)==n[0] || (int)dst.size(0)==n[0]+2*ng[0] );
bool src_halo = (int)src.size(0)==n[0]+2*ng[0];
bool dst_halo = (int)dst.size(0)==n[0]+2*ng[0];
if ( src_halo ) {
ASSERT((int)src.size(0)==nx+2*ngx);
ASSERT((int)src.size(1)==ny+2*ngy);
ASSERT((int)src.size(2)==nz+2*ngz);
ASSERT((int)src.size(0)==n[0]+2*ng[0]);
ASSERT((int)src.size(1)==n[1]+2*ng[1]);
ASSERT((int)src.size(2)==n[2]+2*ng[2]);
} else {
ASSERT((int)src.size(0)==nx);
ASSERT((int)src.size(1)==ny);
ASSERT((int)src.size(2)==nz);
ASSERT((int)src.size(0)==n[0]);
ASSERT((int)src.size(1)==n[1]);
ASSERT((int)src.size(2)==n[2]);
}
if ( dst_halo ) {
ASSERT((int)dst.size(0)==nx+2*ngx);
ASSERT((int)dst.size(1)==ny+2*ngy);
ASSERT((int)dst.size(2)==nz+2*ngz);
ASSERT((int)dst.size(0)==n[0]+2*ng[0]);
ASSERT((int)dst.size(1)==n[1]+2*ng[1]);
ASSERT((int)dst.size(2)==n[2]+2*ng[2]);
} else {
ASSERT((int)dst.size(0)==nx);
ASSERT((int)dst.size(1)==ny);
ASSERT((int)dst.size(2)==nz);
ASSERT((int)dst.size(0)==n[0]);
ASSERT((int)dst.size(1)==n[1]);
ASSERT((int)dst.size(2)==n[2]);
}
if ( src_halo == dst_halo ) {
// Src and dst halos match
@@ -236,19 +287,19 @@ void fillHalo<TYPE>::copy( const Array<TYPE1>& src, Array<TYPE2>& dst )
dst(i) = src(i);
} else if ( src_halo && !dst_halo ) {
// Src has halos
for (int k=0; k<nz; k++) {
for (int j=0; j<ny; j++) {
for (int i=0; i<nx; i++) {
dst(i,j,k) = src(i+ngx,j+ngy,k+ngz);
for (int k=0; k<n[2]; k++) {
for (int j=0; j<n[1]; j++) {
for (int i=0; i<n[0]; i++) {
dst(i,j,k) = src(i+ng[0],j+ng[1],k+ng[2]);
}
}
}
} else if ( !src_halo && dst_halo ) {
// Dst has halos
for (int k=0; k<nz; k++) {
for (int j=0; j<ny; j++) {
for (int i=0; i<nx; i++) {
dst(i+ngx,j+ngy,k+ngz) = src(i,j,k);
for (int k=0; k<n[2]; k++) {
for (int j=0; j<n[1]; j++) {
for (int i=0; i<n[0]; i++) {
dst(i+ng[0],j+ng[1],k+ng[2]) = src(i,j,k);
}
}
}

2
tests/BlobIdentifyParallel.cpp
View File

@@ -100,7 +100,7 @@ int main(int argc, char **argv)
readRankData( rank, nx+2, ny+2, nz+2, Phase, SignDist );
// Communication the halos
fillHalo<double> fillData(comm,rank_info,nx,ny,nz,1,1,1,0,1);
fillHalo<double> fillData(comm,rank_info,{nx,ny,nz},{1,1,1},0,1);
fillData.fill(Phase);
fillData.fill(SignDist);

2
tests/ComponentLabel.cpp
View File

@@ -431,7 +431,7 @@ int main(int argc, char **argv)
Averages.PrintComponents(timestep);
// Create the MeshDataStruct
fillHalo<double> fillData(Dm.Comm,Dm.rank_info,Nx-2,Ny-2,Nz-2,1,1,1,0,1);
fillHalo<double> fillData(Dm.Comm,Dm.rank_info,{Nx-2,Ny-2,Nz-2},{1,1,1},0,1);
std::vector<IO::MeshDataStruct> meshData(1);
meshData[0].meshName = "domain";
meshData[0].mesh = std::shared_ptr<IO::DomainMesh>( new IO::DomainMesh(Dm.rank_info,Nx-2,Ny-2,Nz-2,Lx,Ly,Lz) );

8
tests/TestBlobIdentify.cpp
View File

@@ -133,8 +133,8 @@ void shift_data( DoubleArray& data, int sx, int sy, int sz, const RankInfoStruct
int ngy = ny+2*abs(sy);
int ngz = nz+2*abs(sz);
Array<double> tmp1(nx,ny,nz), tmp2(ngx,ngy,ngz), tmp3(ngx,ngy,ngz);
fillHalo<double> fillData1(comm,rank_info,nx,ny,nz,1,1,1,0,1);
fillHalo<double> fillData2(comm,rank_info,nx,ny,nz,abs(sx),abs(sy),abs(sz),0,1);
fillHalo<double> fillData1(comm,rank_info,{nx,ny,nz},{1,1,1},0,1);
fillHalo<double> fillData2(comm,rank_info,{nx,ny,nz},{abs(sx),abs(sy),abs(sz)},0,1);
fillData1.copy(data,tmp1);
fillData2.copy(tmp1,tmp2);
fillData2.fill(tmp2);
@@ -189,7 +189,7 @@ int main(int argc, char **argv)
fillBubbleData( bubbles, Phase, SignDist, Lx, Ly, Lz, rank_info );
// Communication the halos
fillHalo<double> fillData(comm,rank_info,nx,ny,nz,1,1,1,0,1);
fillHalo<double> fillData(comm,rank_info,{nx,ny,nz},{1,1,1},0,1);
fillData.fill(Phase);
fillData.fill(SignDist);
@@ -205,7 +205,7 @@ int main(int argc, char **argv)
// Create the MeshDataStruct
std::vector<IO::MeshDataStruct> meshData(1);
meshData[0].meshName = "domain";
meshData[0].mesh = std::shared_ptr<IO::DomainMesh>( new IO::DomainMesh(rank_info,nx,ny,nz,Lx,Ly,Lz) );
meshData[0].mesh = std::make_shared<IO::DomainMesh>(rank_info,nx,ny,nz,Lx,Ly,Lz);
std::shared_ptr<IO::Variable> PhaseVar( new IO::Variable() );
std::shared_ptr<IO::Variable> SignDistVar( new IO::Variable() );
std::shared_ptr<IO::Variable> BlobIDVar( new IO::Variable() );

48
tests/TestSegDist.cpp
View File

@@ -11,8 +11,8 @@
#include <fstream>
#include "common/Array.h"
#include "common/Domain.h"
#include "analysis/eikonal.h"
#include "IO/Writer.h"
#include "analysis/distance.h"
std::shared_ptr<Database> loadInputs( int nprocs )
@@ -21,7 +21,7 @@ std::shared_ptr<Database> loadInputs( int nprocs )
auto db = std::make_shared<Database>( );
db->putScalar<int>( "BC", 0 );
db->putVector<int>( "nproc", { 2, 2, 2 } );
db->putVector<int>( "n", { 50, 50, 50 } );
db->putVector<int>( "n", { 100, 100, 100 } );
db->putScalar<int>( "nspheres", 0 );
db->putVector<double>( "L", { 1, 1, 1 } );
return db;
@@ -63,14 +63,13 @@ int main(int argc, char **argv)
int ny = Ny+2;
int nz = Nz+2;
double BubbleRadius = 25;
double BubbleRadius = 0.3*nx;
// Initialize the bubble
double Cx = 1.0*nx;
double Cy = 1.0*ny;
double Cz = 1.0*nz;
DoubleArray Distance(nx,ny,nz);
DoubleArray TrueDist(nx,ny,nz);
Array<char> id(nx,ny,nz);
id.fill(0);
@@ -97,31 +96,15 @@ int main(int argc, char **argv)
}
}
// Initialize the signed distance function
for (int k=0;k<nz;k++){
for (int j=0;j<ny;j++){
for (int i=0;i<nx;i++){
// Initialize distance to +/- 1
Distance(i,j,k) = 2.0*id(i,j,k)-1.0;
}
}
}
if (rank==0) printf("Nx = %i \n",(int)Distance.size(0));
if (rank==0) printf("Ny = %i \n",(int)Distance.size(1));
if (rank==0) printf("Nz = %i \n",(int)Distance.size(2));
MPI_Barrier(comm);
if (rank==0) printf("Initialized! Converting to Signed Distance function \n");
double starttime,stoptime,cputime;
starttime = MPI_Wtime();
double err1 = Eikonal(Distance,id,Dm,1000);
stoptime = MPI_Wtime();
cputime = (stoptime - starttime);
if (rank==0) printf("Total time: %f seconds \n",cputime);
double t1 = MPI_Wtime();
DoubleArray Distance(nx,ny,nz);
CalcDist(Distance,id,Dm,{false,false,false});
double t2 = MPI_Wtime();
if (rank==0)
printf("Total time: %f seconds \n",t2-t1);
double localError=0.0;
int localCount = 0;
@@ -140,7 +123,6 @@ int main(int argc, char **argv)
MPI_Allreduce(&localError,&globalError,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
MPI_Allreduce(&localCount,&globalCount,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
double err2 = sqrt(globalError)/(double (globalCount));
if (rank==0) printf("global variation %f \n", err1);
if (rank==0) printf("Mean error %f \n", err2);
// Write the results to visit
@@ -149,17 +131,17 @@ int main(int argc, char **argv)
Array<double> ID(Nx,Ny,Nz);
Array<double> dist1(Nx,Ny,Nz);
Array<double> dist2(Nx,Ny,Nz);
fillHalo<double> fillData(Dm.Comm, Dm.rank_info,Nx,Ny,Nz,1,1,1,0,1);
fillHalo<double> fillData(Dm.Comm, Dm.rank_info,{Nx,Ny,Nz},{1,1,1},0,1);
fillData.copy( ID0, ID );
fillData.copy( Distance, dist1 );
fillData.copy( TrueDist, dist2 );
fillData.copy( TrueDist, dist1 );
fillData.copy( Distance, dist2 );
std::vector<IO::MeshDataStruct> data(1);
data[0].meshName = "mesh";
data[0].mesh.reset( new IO::DomainMesh( Dm.rank_info, Nx, Ny, Nz, Dm.Lx, Dm.Ly, Dm.Lz ) );
data[0].vars.emplace_back( new IO::Variable( 1, IO::VariableType::VolumeVariable, "ID", ID ) );
data[0].vars.emplace_back( new IO::Variable( 1, IO::VariableType::VolumeVariable, "Distance", dist1 ) );
data[0].vars.emplace_back( new IO::Variable( 1, IO::VariableType::VolumeVariable, "TrueDist", dist2 ) );
data[0].vars.emplace_back( new IO::Variable( 1, IO::VariableType::VolumeVariable, "error", dist1-dist2 ) );
data[0].vars.emplace_back( new IO::Variable( 1, IO::VariableType::VolumeVariable, "TrueDist", dist1 ) );
data[0].vars.emplace_back( new IO::Variable( 1, IO::VariableType::VolumeVariable, "Distance", dist2 ) );
data[0].vars.emplace_back( new IO::Variable( 1, IO::VariableType::VolumeVariable, "error", dist2-dist1 ) );
IO::initialize( "", "silo", false );
IO::writeData( "testSegDist", data, MPI_COMM_WORLD );

72
tests/lbpm_morph_pp.cpp
View File

@@ -343,42 +343,42 @@ int main(int argc, char **argv)
PackID(Dm.sendList_yZ, Dm.sendCount_yZ ,sendID_yZ, id);
PackID(Dm.sendList_YZ, Dm.sendCount_YZ ,sendID_YZ, id);
//......................................................................................
MPI_Sendrecv(sendID_x,Dm.sendCount_x,MPI_CHAR,Dm.rank_x,sendtag,
recvID_X,Dm.recvCount_X,MPI_CHAR,Dm.rank_X,recvtag,comm,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_X,Dm.sendCount_X,MPI_CHAR,Dm.rank_X,sendtag,
recvID_x,Dm.recvCount_x,MPI_CHAR,Dm.rank_x,recvtag,comm,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_y,Dm.sendCount_y,MPI_CHAR,Dm.rank_y,sendtag,
recvID_Y,Dm.recvCount_Y,MPI_CHAR,Dm.rank_Y,recvtag,comm,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_Y,Dm.sendCount_Y,MPI_CHAR,Dm.rank_Y,sendtag,
recvID_y,Dm.recvCount_y,MPI_CHAR,Dm.rank_y,recvtag,comm,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_z,Dm.sendCount_z,MPI_CHAR,Dm.rank_z,sendtag,
recvID_Z,Dm.recvCount_Z,MPI_CHAR,Dm.rank_Z,recvtag,comm,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_Z,Dm.sendCount_Z,MPI_CHAR,Dm.rank_Z,sendtag,
recvID_z,Dm.recvCount_z,MPI_CHAR,Dm.rank_z,recvtag,comm,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_xy,Dm.sendCount_xy,MPI_CHAR,Dm.rank_xy,sendtag,
recvID_XY,Dm.recvCount_XY,MPI_CHAR,Dm.rank_XY,recvtag,comm,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_XY,Dm.sendCount_XY,MPI_CHAR,Dm.rank_XY,sendtag,
recvID_xy,Dm.recvCount_xy,MPI_CHAR,Dm.rank_xy,recvtag,comm,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_Xy,Dm.sendCount_Xy,MPI_CHAR,Dm.rank_Xy,sendtag,
recvID_xY,Dm.recvCount_xY,MPI_CHAR,Dm.rank_xY,recvtag,comm,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_xY,Dm.sendCount_xY,MPI_CHAR,Dm.rank_xY,sendtag,
recvID_Xy,Dm.recvCount_Xy,MPI_CHAR,Dm.rank_Xy,recvtag,comm,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_xz,Dm.sendCount_xz,MPI_CHAR,Dm.rank_xz,sendtag,
recvID_XZ,Dm.recvCount_XZ,MPI_CHAR,Dm.rank_XZ,recvtag,comm,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_XZ,Dm.sendCount_XZ,MPI_CHAR,Dm.rank_XZ,sendtag,
recvID_xz,Dm.recvCount_xz,MPI_CHAR,Dm.rank_xz,recvtag,comm,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_Xz,Dm.sendCount_Xz,MPI_CHAR,Dm.rank_Xz,sendtag,
recvID_xZ,Dm.recvCount_xZ,MPI_CHAR,Dm.rank_xZ,recvtag,comm,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_xZ,Dm.sendCount_xZ,MPI_CHAR,Dm.rank_xZ,sendtag,
recvID_Xz,Dm.recvCount_Xz,MPI_CHAR,Dm.rank_Xz,recvtag,comm,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_yz,Dm.sendCount_yz,MPI_CHAR,Dm.rank_yz,sendtag,
recvID_YZ,Dm.recvCount_YZ,MPI_CHAR,Dm.rank_YZ,recvtag,comm,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_YZ,Dm.sendCount_YZ,MPI_CHAR,Dm.rank_YZ,sendtag,
recvID_yz,Dm.recvCount_yz,MPI_CHAR,Dm.rank_yz,recvtag,comm,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_Yz,Dm.sendCount_Yz,MPI_CHAR,Dm.rank_Yz,sendtag,
recvID_yZ,Dm.recvCount_yZ,MPI_CHAR,Dm.rank_yZ,recvtag,comm,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_yZ,Dm.sendCount_yZ,MPI_CHAR,Dm.rank_yZ,sendtag,
recvID_Yz,Dm.recvCount_Yz,MPI_CHAR,Dm.rank_Yz,recvtag,comm,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_x,Dm.sendCount_x,MPI_CHAR,Dm.rank_x(),sendtag,
recvID_X,Dm.recvCount_X,MPI_CHAR,Dm.rank_X(),recvtag,comm,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_X,Dm.sendCount_X,MPI_CHAR,Dm.rank_X(),sendtag,
recvID_x,Dm.recvCount_x,MPI_CHAR,Dm.rank_x(),recvtag,comm,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_y,Dm.sendCount_y,MPI_CHAR,Dm.rank_y(),sendtag,
recvID_Y,Dm.recvCount_Y,MPI_CHAR,Dm.rank_Y(),recvtag,comm,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_Y,Dm.sendCount_Y,MPI_CHAR,Dm.rank_Y(),sendtag,
recvID_y,Dm.recvCount_y,MPI_CHAR,Dm.rank_y(),recvtag,comm,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_z,Dm.sendCount_z,MPI_CHAR,Dm.rank_z(),sendtag,
recvID_Z,Dm.recvCount_Z,MPI_CHAR,Dm.rank_Z(),recvtag,comm,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_Z,Dm.sendCount_Z,MPI_CHAR,Dm.rank_Z(),sendtag,
recvID_z,Dm.recvCount_z,MPI_CHAR,Dm.rank_z(),recvtag,comm,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_xy,Dm.sendCount_xy,MPI_CHAR,Dm.rank_xy(),sendtag,
recvID_XY,Dm.recvCount_XY,MPI_CHAR,Dm.rank_XY(),recvtag,comm,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_XY,Dm.sendCount_XY,MPI_CHAR,Dm.rank_XY(),sendtag,
recvID_xy,Dm.recvCount_xy,MPI_CHAR,Dm.rank_xy(),recvtag,comm,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_Xy,Dm.sendCount_Xy,MPI_CHAR,Dm.rank_Xy(),sendtag,
recvID_xY,Dm.recvCount_xY,MPI_CHAR,Dm.rank_xY(),recvtag,comm,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_xY,Dm.sendCount_xY,MPI_CHAR,Dm.rank_xY(),sendtag,
recvID_Xy,Dm.recvCount_Xy,MPI_CHAR,Dm.rank_Xy(),recvtag,comm,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_xz,Dm.sendCount_xz,MPI_CHAR,Dm.rank_xz(),sendtag,
recvID_XZ,Dm.recvCount_XZ,MPI_CHAR,Dm.rank_XZ(),recvtag,comm,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_XZ,Dm.sendCount_XZ,MPI_CHAR,Dm.rank_XZ(),sendtag,
recvID_xz,Dm.recvCount_xz,MPI_CHAR,Dm.rank_xz(),recvtag,comm,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_Xz,Dm.sendCount_Xz,MPI_CHAR,Dm.rank_Xz(),sendtag,
recvID_xZ,Dm.recvCount_xZ,MPI_CHAR,Dm.rank_xZ(),recvtag,comm,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_xZ,Dm.sendCount_xZ,MPI_CHAR,Dm.rank_xZ(),sendtag,
recvID_Xz,Dm.recvCount_Xz,MPI_CHAR,Dm.rank_Xz(),recvtag,comm,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_yz,Dm.sendCount_yz,MPI_CHAR,Dm.rank_yz(),sendtag,
recvID_YZ,Dm.recvCount_YZ,MPI_CHAR,Dm.rank_YZ(),recvtag,comm,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_YZ,Dm.sendCount_YZ,MPI_CHAR,Dm.rank_YZ(),sendtag,
recvID_yz,Dm.recvCount_yz,MPI_CHAR,Dm.rank_yz(),recvtag,comm,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_Yz,Dm.sendCount_Yz,MPI_CHAR,Dm.rank_Yz(),sendtag,
recvID_yZ,Dm.recvCount_yZ,MPI_CHAR,Dm.rank_yZ(),recvtag,comm,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_yZ,Dm.sendCount_yZ,MPI_CHAR,Dm.rank_yZ(),sendtag,
recvID_Yz,Dm.recvCount_Yz,MPI_CHAR,Dm.rank_Yz(),recvtag,comm,MPI_STATUS_IGNORE);
UnpackID(Dm.recvList_x, Dm.recvCount_x ,recvID_x, id);
UnpackID(Dm.recvList_X, Dm.recvCount_X ,recvID_X, id);

2
tests/lbpm_morphdrain_pp.cpp
View File

@@ -155,7 +155,7 @@ int main(int argc, char **argv)
xdim=Dm.Nx-2;
ydim=Dm.Ny-2;
zdim=Dm.Nz-2;
fillHalo<double> fillData(Dm.Comm, Dm.rank_info,xdim,ydim,zdim,1,1,1,0,1);
fillHalo<double> fillData(Dm.Comm, Dm.rank_info,{xdim,ydim,zdim},{1,1,1},0,1);
DoubleArray SignDist(nx,ny,nz);
// Read the signed distance from file

2
tests/lbpm_refine_pp.cpp
View File

@@ -87,7 +87,7 @@ int main(int argc, char **argv)
// Communication the halos
const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz);
fillHalo<double> fillData(comm,rank_info,rnx,rny,rnz,1,1,1,0,1);
fillHalo<double> fillData(comm,rank_info,{rnx,rny,rnz},{1,1,1},0,1);
nx+=2; ny+=2; nz+=2;
rnx+=2; rny+=2; rnz+=2;

11
tests/lbpm_segmented_pp.cpp
View File

@@ -12,7 +12,7 @@
#include "common/Array.h"
#include "common/Domain.h"
#include "analysis/TwoPhase.h"
#include "analysis/eikonal.h"
#include "analysis/distance.h"
inline void MeanFilter(DoubleArray &Mesh){
for (int k=1; k<(int)Mesh.size(2)-1; k++){
@@ -251,15 +251,8 @@ int main(int argc, char **argv)
}
MeanFilter(Averages.SDs);
double LocalVar, TotalVar;
if (rank==0) printf("Initialized solid phase -- Converting to Signed Distance function \n");
int Maxtime=10*max(max(Dm.Nx*Dm.nprocx(),Dm.Ny*Dm.nprocy()),Dm.Nz*Dm.nprocz());
Maxtime=min(Maxtime,MAXTIME);
LocalVar = Eikonal(Averages.SDs,id,Dm,Maxtime);
MPI_Allreduce(&LocalVar,&TotalVar,1,MPI_DOUBLE,MPI_SUM,comm);
TotalVar /= nprocs;
if (rank==0) printf("Final variation in signed distance function %f \n",TotalVar);
CalcDist(Averages.SDs,id,Dm);
sprintf(LocalRankFilename,"SignDist.%05i",rank);
FILE *DIST = fopen(LocalRankFilename,"wb");

2
tests/testCommunication.cpp
View File

@@ -221,7 +221,7 @@ int testHalo( MPI_Comm comm, int nprocx, int nprocy, int nprocz, int depth )
}
// Communicate the halo
fillHalo<TYPE> fillData(comm,rank_info,nx,ny,nz,1,1,1,0,depth);
fillHalo<TYPE> fillData(comm,rank_info,{nx,ny,nz},{1,1,1},0,depth);
fillData.fill(array);
// Check the results