added subphase test
This commit is contained in:
JamesEMcclure
parent
07f4356869
commit
b00c382206
146
tests/TestSubphase.cpp
Normal file
146
tests/TestSubphase.cpp
Normal file
@ -0,0 +1,146 @@
|
||||
// Sequential blob analysis
|
||||
// Reads parallel simulation data and performs connectivity analysis
|
||||
// and averaging on a blob-by-blob basis
|
||||
// James E. McClure 2014
|
||||
|
||||
#include <iostream>
|
||||
#include <math.h>
|
||||
#include "common/Communication.h"
|
||||
#include "analysis/analysis.h"
|
||||
#include "analysis/SubPhase.h"
|
||||
|
||||
|
||||
std::shared_ptr<Database> loadInputs( int nprocs )
|
||||
{
|
||||
//auto db = std::make_shared<Database>( "Domain.in" );
|
||||
auto db = std::make_shared<Database>();
|
||||
db->putScalar<int>( "BC", 0 );
|
||||
db->putVector<int>( "nproc", { 1, 1, 1 } );
|
||||
db->putVector<int>( "n", { 100, 100, 100 } );
|
||||
db->putScalar<int>( "nspheres", 1 );
|
||||
db->putVector<double>( "L", { 1, 1, 1 } );
|
||||
return db;
|
||||
}
|
||||
|
||||
|
||||
int main(int argc, char **argv)
|
||||
{
|
||||
// Initialize MPI
|
||||
int rank, nprocs;
|
||||
MPI_Init(&argc,&argv);
|
||||
MPI_Comm comm = MPI_COMM_WORLD;
|
||||
MPI_Comm_rank(comm,&rank);
|
||||
MPI_Comm_size(comm,&nprocs);
|
||||
{ // Limit scope so variables that contain communicators will free before MPI_Finialize
|
||||
|
||||
if ( rank==0 ) {
|
||||
printf("-----------------------------------------------------------\n");
|
||||
printf("Unit test for torus (Euler-Poincarie characteristic) \n");
|
||||
printf("-----------------------------------------------------------\n");
|
||||
}
|
||||
|
||||
//.......................................................................
|
||||
// Reading the domain information file
|
||||
//.......................................................................
|
||||
int i,j,k,n;
|
||||
|
||||
// Load inputs
|
||||
auto db = loadInputs( nprocs );
|
||||
int Nx = db->getVector<int>( "n" )[0];
|
||||
int Ny = db->getVector<int>( "n" )[1];
|
||||
int Nz = db->getVector<int>( "n" )[2];
|
||||
int nprocx = db->getVector<int>( "nproc" )[0];
|
||||
int nprocy = db->getVector<int>( "nproc" )[1];
|
||||
int nprocz = db->getVector<int>( "nproc" )[2];
|
||||
|
||||
if (rank==0){
|
||||
printf("********************************************************\n");
|
||||
printf("Sub-domain size = %i x %i x %i\n",Nx,Ny,Nz);
|
||||
printf("********************************************************\n");
|
||||
}
|
||||
|
||||
// Get the rank info
|
||||
std::shared_ptr<Domain> Dm(new Domain(db,comm));
|
||||
// const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz);
|
||||
std::shared_ptr<SubPhase> Averages(new SubPhase(Dm));
|
||||
Nx += 2;
|
||||
Ny += 2;
|
||||
Nz += 2;
|
||||
int N = Nx*Ny*Nz;
|
||||
//.......................................................................
|
||||
for ( k=1;k<Nz-1;k++){
|
||||
for ( j=1;j<Ny-1;j++){
|
||||
for ( i=1;i<Nx-1;i++){
|
||||
n = k*Nx*Ny+j*Nx+i;
|
||||
Dm->id[n] = 1;
|
||||
}
|
||||
}
|
||||
}
|
||||
//.......................................................................
|
||||
Dm->CommInit(); // Initialize communications for domains
|
||||
//.......................................................................
|
||||
|
||||
//.......................................................................
|
||||
// Assign the phase ID field based and the signed distance
|
||||
//.......................................................................
|
||||
double R1,R2;
|
||||
double CX,CY,CZ; //CY1,CY2;
|
||||
CX=Nx*nprocx*0.5;
|
||||
CY=Ny*nprocy*0.5;
|
||||
CZ=Nz*nprocz*0.5;
|
||||
R1 = Nx*nprocx*0.2; // middle radius
|
||||
R2 = Nx*nprocx*0.1; // donut thickness
|
||||
//
|
||||
//CY1=Nx*nprocx*0.5+R1;
|
||||
//CY2=Ny*nprocy*0.5-R1;
|
||||
|
||||
double x,y,z;
|
||||
if (rank==0) printf("Initializing the system \n");
|
||||
for ( k=1;k<Nz-1;k++){
|
||||
for ( j=1;j<Ny-1;j++){
|
||||
for ( i=1;i<Nx-1;i++){
|
||||
n = k*Nx*Ny+j*Nx+i;
|
||||
|
||||
// global position relative to center
|
||||
x = Dm->iproc()*Nx+i - CX;
|
||||
y = Dm->jproc()*Ny+j - CY;
|
||||
z = Dm->kproc()*Nz+k - CZ;
|
||||
|
||||
// Shrink the sphere sizes by two voxels to make sure they don't touch
|
||||
Averages->SDs(i,j,k) = 100.0;
|
||||
//..............................................................................
|
||||
// Single torus
|
||||
Averages->Phi(i,j,k) = sqrt((sqrt(x*x+y*y) - R1)*(sqrt(x*x+y*y) - R1) + z*z) - R2;
|
||||
// Double torus
|
||||
/* y = Dm->jproc()*Ny+j - CY1;
|
||||
//z = Dm->kproc()*Nz+k - CZ +R1;
|
||||
Averages->Phi(i,j,k) = sqrt((sqrt(x*x+y*y) - R1)*(sqrt(x*x+y*y) - R1) + z*z) - R2;
|
||||
|
||||
y = Dm->jproc()*Ny+j - CY2;
|
||||
//z = Dm->kproc()*Nz+k - CZ-R1;
|
||||
Averages->Phi(i,j,k) = min(Averages->Phi(i,j,k),
|
||||
sqrt((sqrt(x*x+y*y) - R1)*(sqrt(x*x+y*y) - R1) + z*z) - R2);
|
||||
*///..............................................................................
|
||||
|
||||
//Averages->Phi(i,j,k) = - Averages->Phi(i,j,k);
|
||||
if (Averages->Phi(i,j,k) > 0.0){
|
||||
Dm->id[n] = 2;
|
||||
}
|
||||
else{
|
||||
Dm->id[n] = 1;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (rank==0) printf("Aggregating labels \n");
|
||||
|
||||
Averages->AggregateLabels("phase_id.raw");
|
||||
// Averages->Reduce();
|
||||
|
||||
} // Limit scope so variables that contain communicators will free before MPI_Finialize
|
||||
MPI_Barrier(comm);
|
||||
MPI_Finalize();
|
||||
return 0;
|
||||
}
|
||||
|
||||
Loading…
Reference in New Issue
Block a user