/****************************************************************/
/* */
/* Program sinter3d.c */
/* To simulate 3-D curvature controlled sintering */
/* of a loosely packed powder */
/* This version to sinter between two phases in */
/* a multi-phase system */
/* This version is periodic in all three directions */
/* Application: Modify local structure of generated */
/* 3-D porous media by matching a specific hydraulic */
/* radius */
/* */
/* Programmer: Dale P. Bentz */
/* NIST */
/* 100 Bureau Drive Mail Stop 8621 */
/* Gaithersburg, MD 20899-8621 */
/* Phone: (301) 975-5865 */
/* E-mail: dale.bentz@.nist.gov */
/* Date: Summer 1994 */
/* */
/****************************************************************/
#include <stdio.h>
#include <math.h>
#define SYSSIZE 100
#define MAXCYC 1000 /* maximum sintering cycles to use */
#define MAXSPH 10000 /* maximum number of elements in a spherical template */
/* array phase stores the microstructure representation */
/* array curvature stores the local curvature values */
static unsigned short int phase [SYSSIZE+1] [SYSSIZE+1] [SYSSIZE+1];
static unsigned short int curvature [SYSSIZE+1] [SYSSIZE+1] [SYSSIZE+1];
int nsph,xsph[MAXSPH],ysph[MAXSPH],zsph[MAXSPH];
int *seed;
long int nsolid[500],nair[500];
/* Random number generator ran1 from Computers in Physics */
/* Volume 6 No. 5, 1992, 522-524, Press and Teukolsky */
/* To generate real random numbers 0.0-1.0 */
/* Should be seeded with a negative integer */
#define IA 16807
#define IM 2147483647
#define IQ 127773
#define IR 2836
#define NTAB 32
#define EPS (1.2E-07)
#define MAX(a,b) (a>b)?a:b
#define MIN(a,b) (a<b)?a:b
double ran1(idum)
int *idum;
/* Calls: no routines */
/* Called by: gsphere,makefloc */
{
int j,k;
static int iv[NTAB],iy=0;
void nrerror();
static double NDIV = 1.0/(1.0+(IM-1.0)/NTAB);
static double RNMX = (1.0-EPS);
static double AM = (1.0/IM);
if ((*idum <= 0) || (iy == 0)) {
*idum = MAX(-*idum,*idum);
for(j=NTAB+7;j>=0;j--) {
k = *idum/IQ;
*idum = IA*(*idum-k*IQ)-IR*k;
if(*idum < 0) *idum += IM;
if(j < NTAB) iv[j] = *idum;
}
iy = iv[0];
}
k = *idum/IQ;
*idum = IA*(*idum-k*IQ)-IR*k;
if(*idum<0) *idum += IM;
j = iy*NDIV;
iy = iv[j];
iv[j] = *idum;
return MIN(AM*iy,RNMX);
}
#undef IA
#undef IM
#undef IQ
#undef IR
#undef NTAB
#undef EPS
#undef MAX
#undef MIN
/* routine to create a template for the sphere of interest of radius size */
/* to be used in curvature evaluation */
/* Called by: runsint */
/* Calls no other routines */
int maketemp(size)
int size;
{
int icirc,xval,yval,zval;
float xtmp,ytmp;
float dist;
/* determine and store the locations of all pixels in the 3-D sphere */
icirc=0;
for(xval=(-size);xval<=size;xval++){
xtmp=(float)(xval*xval);
for(yval=(-size);yval<=size;yval++){
ytmp=(float)(yval*yval);
for(zval=(-size);zval<=size;zval++){
dist=sqrt(xtmp+ytmp+(float)(zval*zval));
if(dist<=((float)size+0.5)){
icirc+=1;
if(icirc>=MAXSPH){
printf("Too many elements in sphere \n");
printf("Must change value of MAXSPH parameter \n");
printf("Currently set at %d \n",MAXSPH);
exit(1);
}
xsph[icirc]=xval;
ysph[icirc]=yval;
zsph[icirc]=zval;
}
}
}
}
/* return the number of pixels contained in sphere of radius (size+0.5) */
return(icirc);
}
/* routine to count phase fractions (porosity and solids) */
/* Called by main routine */
/* Calls no other routines */
void phcount()
{
long int npore,nsolid [37];
int ix,iy,iz;
npore=0;
for(ix=1;ix<37;ix++){
nsolid[ix]=0;
}
/* check all pixels in the 3-D system */
for(ix=0;ix<SYSSIZE;ix++){
for(iy=0;iy<SYSSIZE;iy++){
for(iz=0;iz<SYSSIZE;iz++){
if(phase [ix] [iy] [iz]==0){
npore+=1;
}
else{
nsolid[phase [ix] [iy] [iz]]+=1;
}
}
}
}
printf("Pores are: %ld \n",npore);
printf("Solids are: %ld %ld %ld %ld %ld %ld\n",nsolid[1],nsolid[2],
nsolid[3],nsolid[4],nsolid[5],nsolid[6]);
}
/* routine to return number of surface faces exposed to porosity */
/* for pixel located at (xin,yin,zin) */
/* Called by rhcalc */
/* Calls no other routines */
int surfpix(xin,yin,zin)
int xin,yin,zin;
{
int npix,ix1,iy1,iz1;
npix=0;
/* check each of the six immediate neighbors */
/* using periodic boundary conditions */
ix1=xin-1;
if(ix1<0){ix1+=SYSSIZE;}
if(phase[ix1][yin][zin]==0){
npix+=1;
}
ix1=xin+1;
if(ix1>=SYSSIZE){ix1-=SYSSIZE;}
if(phase[ix1][yin][zin]==0){
npix+=1;
}
iy1=yin-1;
if(iy1<0){iy1+=SYSSIZE;}
if(phase[xin][iy1][zin]==0){
npix+=1;
}
iy1=yin+1;
if(iy1>=SYSSIZE){iy1-=SYSSIZE;}
if(phase[xin][iy1][zin]==0){
npix+=1;
}
iz1=zin-1;
if(iz1<0){iz1+=SYSSIZE;}
if(phase[xin][yin][iz1]==0){
npix+=1;
}
iz1=zin+1;
if(iz1>=SYSSIZE){iz1-=SYSSIZE;}
if(phase[xin][yin][iz1]==0){
npix+=1;
}
return(npix);
}
/* routine to return the current hydraulic radius for phase phin */
/* Calls surfpix */
/* Called by runsint */
float rhcalc(phin)
int phin;
{
int ix,iy,iz;
long int porc,surfc;
float rhval;
porc=surfc=0;
/* Check all pixels in the 3-D volume */
for(ix=0;ix<SYSSIZE;ix++){
for(iy=0;iy<SYSSIZE;iy++){
for(iz=0;iz<SYSSIZE;iz++){
if(phase [ix] [iy] [iz]==phin){
porc+=1;
surfc+=surfpix(ix,iy,iz);
}
}
}
}
printf("Phase area count is %ld \n",porc);
printf("Phase surface count is %ld \n",surfc);
rhval=(float)porc*6./(4.*(float)surfc);
printf("Hydraulic radius is %f \n",rhval);
return(rhval);
}
/* routine to return count of pixels in a spherical template which are phase */
/* phin or porosity (phase=0) */
/* Calls no other routines */
/* Called by sysinit */
int countem(xp,yp,zp,phin)
int xp,yp,zp,phin;
{
int xc,yc,zc;
int cumnum,ic;
cumnum=0;
for(ic=1;ic<=nsph;ic++){
xc=xp+xsph[ic];
yc=yp+ysph[ic];
zc=zp+zsph[ic];
/* Use periodic boundaries */
if(xc<0){xc+=SYSSIZE;}
else if(xc>=SYSSIZE){xc-=SYSSIZE;}
if(yc<0){yc+=SYSSIZE;}
else if(yc>=SYSSIZE){yc-=SYSSIZE;}
if(zc<0){zc+=SYSSIZE;}
else if(zc>=SYSSIZE){zc-=SYSSIZE;}
if((xc!=xp)||(yc!=yp)||(zc!=zp)){
if((phase [xc] [yc] [zc]==phin)||(phase [xc] [yc] [zc]==0)){
cumnum+=1;
}
}
}
return(cumnum);
}
/* routine to initialize system by determining local curvature */
/* of all phase 1 and phase 2 pixels */
/* Calls countem */
/* Called by runsint */
void sysinit(ph1,ph2)
int ph1,ph2;
{
int count,xl,yl,zl;
count=0;
/* process all pixels in the 3-D box */
for(xl=0;xl<SYSSIZE;xl++){
for(yl=0;yl<SYSSIZE;yl++){
for(zl=0;zl<SYSSIZE;zl++){
/* determine local curvature */
/* For phase 1 want to determine number of porosity pixels */
/* (phase=0) in immediate neighborhood */
if(phase [xl] [yl] [zl]==ph1){
count=countem(xl,yl,zl,0);
}
/* For phase 2 want to determine number of porosity or phase */
/* 2 pixels in immediate neighborhood */
if(phase [xl] [yl] [zl]==ph2){
count=countem(xl,yl,zl,ph2);
}
if((count<0)||(count>=nsph)){
printf("Error count is %d \n",count);
printf("xl %d yl %d zl %d \n",xl,yl,zl);
}
/* case where we have a phase 1 surface pixel */
/* with non-zero local curvature */
if((count>=0)&&(phase [xl] [yl] [zl]==ph1)){
curvature [xl] [yl] [zl]=count;
/* update solid curvature histogram */
nsolid[count]+=1;
}
/* case where we have a phase 2 surface pixel */
if((count>=0)&&(phase [xl] [yl] [zl]==ph2)){
curvature [xl] [yl] [zl]=count;
/* update air curvature histogram */
nair[count]+=1;
}
}
}
} /* end of xl loop */
}
/* routine to scan system and determine nsolid (ph2) and nair (ph1) */
/* histograms based on values in phase and curvature arrays */
/* Calls no other routines */
/* Called by runsint */
void sysscan(ph1,ph2)
int ph1,ph2;
{
int xd,yd,zd,curvval;
/* Scan all pixels in 3-D system */
for(xd=0;xd<SYSSIZE;xd++){
for(yd=0;yd<SYSSIZE;yd++){
for(zd=0;zd<SYSSIZE;zd++){
curvval=curvature [xd] [yd] [zd];
if(phase [xd] [yd] [zd]==ph2){
nair[curvval]+=1;
}
else if (phase [xd] [yd] [zd]==ph1){
nsolid[curvval]+=1;
}
}
}
}
}
/* routine to return how many cells of solid curvature histogram to use */
/* to accomodate nsearch pixels moving */
/* want to use highest values first */
/* Calls no other routines */
/* Called by movepix */
int procsol(nsearch)
int nsearch;
{
int valfound,i,stop;
long int nsofar;
/* search histogram from top down until cumulative count */
/* exceeds nsearch */
valfound=nsph-1;
nsofar=0;
stop=0;
for(i=(nsph-1);((i>=0)&&(stop==0));i--){
nsofar+=nsolid[i];
if(nsofar>nsearch){
valfound=i;
stop=1;
}
}
return(valfound);
}
/* routine to determine how many cells of air curvature histogram to use */
/* to accomodate nsearch moving pixels */
/* want to use lowest values first */
/* Calls no other routines */
/* Called by movepix */
int procair(nsearch)
int nsearch;
{
int valfound,i,stop;
long int nsofar;
/* search histogram from bottom up until cumulative count */
/* exceeds nsearch */
valfound=0;
nsofar=0;
stop=0;
for(i=0;((i<nsph)&&(stop==0));i++){
nsofar+=nair[i];
if(nsofar>nsearch){
valfound=i;
stop=1;
}
}
return(valfound);
}
/* routine to move requested number of pixels (ntomove) from highest */
/* curvature phase 1 (ph1) sites to lowest curvature phase 2 (ph2) sites */
/* Calls procsol and procair */
/* Called by runsint */
int movepix(ntomove,ph1,ph2)
int ntomove,ph1,ph2;
{
int xloc[2100],yloc[2100],zloc[2100];
int count1,count2,ntot,countc,i,xp,yp,zp;
int cmin,cmax,cfg;
int alldone;
long int nsolc,nairc,nsum,nsolm,nairm,nst1,nst2,next1,next2;
float pck,plsol,plair;
alldone=0;
/* determine critical values for removal and placement */
count1=procsol(ntomove);
nsum=0;
cfg=0;
cmax=count1;
for(i=nsph;i>count1;i--){
if((nsolid[i]>0)&&(cfg==0)){
cfg=1;
cmax=i;
}
nsum+=nsolid[i];
}
/* Determine movement probability for last cell */
plsol=(float)(ntomove-nsum)/(float)nsolid[count1];
next1=ntomove-nsum;
nst1=nsolid[count1];
count2=procair(ntomove);
nsum=0;
cmin=count2;
cfg=0;
for(i=0;i<count2;i++){
if((nair[i]>0)&&(cfg==0)){
cfg=1;
cmin=i;
}
nsum+=nair[i];
}
/* Determine movement probability for last cell */
plair=(float)(ntomove-nsum)/(float)nair[count2];
next2=ntomove-nsum;
nst2=nair[count2];
/* Check to see if equilibrium has been reached --- */
/* no further increase in hydraulic radius is possible */
if(cmin>=cmax){
alldone=1;
printf("Stopping - at equilibrium \n");
printf("cmin- %d cmax- %d \n",cmin,cmax);
return(alldone);
}
/* initialize counters for performing sintering */
ntot=0;
nsolc=0;
nairc=0;
nsolm=0;
nairm=0;
/* Now process each pixel in turn */
for(xp=0;xp<SYSSIZE;xp++){
for(yp=0;yp<SYSSIZE;yp++){
for(zp=0;zp<SYSSIZE;zp++){
countc=curvature [xp] [yp] [zp];
/* handle phase 1 case first */
if(phase [xp] [yp] [zp]==ph1){
if(countc>count1){
/* convert from phase 1 to phase 2 */
phase [xp] [yp] [zp]=ph2;
/* update appropriate histogram cells */
nsolid[countc]-=1;
nair[countc]+=1;
/* store the location of the modified pixel */
ntot+=1;
xloc[ntot]=xp;
yloc[ntot]=yp;
zloc[ntot]=zp;
}
if(countc==count1){
nsolm+=1;
/* generate probability for pixel being removed */
pck=ran1(seed);
if((pck<0)||(pck>1.0)){pck=1.0;}
if(((pck<plsol)&&(nsolc<next1))||((nst1-nsolm)<(next1-nsolc))){
nsolc+=1;
/* convert phase 1 pixel to phase 2 */
phase [xp] [yp] [zp]=ph2;
/* update appropriate histogram cells */
nsolid[count1]-=1;
nair[count1]+=1;
/* store the location of the modified pixel */
ntot+=1;
xloc[ntot]=xp;
yloc[ntot]=yp;
zloc[ntot]=zp;
}
}
}
/* handle phase 2 case here */
else if (phase [xp] [yp] [zp]==ph2){
if(countc<count2){
/* convert phase 2 pixel to phase 1 */
phase [xp] [yp] [zp]=ph1;
nsolid[countc]+=1;
nair[countc]-=1;
ntot+=1;
xloc[ntot]=xp;
yloc[ntot]=yp;
zloc[ntot]=zp;
}
if(countc==count2){
nairm+=1;
pck=ran1(seed);
if((pck<0)||(pck>1.0)){pck=1.0;}
if(((pck<plair)&&(nairc<next2))||((nst2-nairm)<(next2-nairc))){
nairc+=1;
/* convert phase 2 to phase 1 */
phase [xp] [yp] [zp]=ph1;
nsolid[count2]+=1;
nair[count2]-=1;
ntot+=1;
xloc[ntot]=xp;
yloc[ntot]=yp;
zloc[ntot]=zp;
}
}
}
} /* end of zp loop */
} /* end of yp loop */
} /* end of xloop */
printf("ntot is %d \n",ntot);
return(alldone);
}
/* routine to execute user input number of cycles of sintering algorithm */
/* Calls maketemp, rhcalc, sysinit, sysscan, and movepix */
/* Called by main routine */
void runsint()
{
int natonce,ncyc,i,rade,j,rflag;
int ph1id,ph2id,keepgo;
long int curvsum1,curvsum2,pixsum1,pixsum2;
float rhnow,rhtarget,avecurv1,avecurv2;
/* initialize the solid and air count histograms */
for(i=0;i<=499;i++){
nsolid[i]=0;
nair[i]=0;
}
/* Obtain needed user input */
printf("Enter phases to execute sintering between \n");
scanf("%d %d",&ph1id,&ph2id);
printf("%d %d \n",ph1id,ph2id);
printf("Enter number of pixels to move at once (e.g. 200)\n");
scanf("%d",&natonce);
printf("%d \n",natonce);
printf("Enter target hydraulic radius \n");
scanf("%f",&rhtarget);
printf("%f \n",rhtarget);
printf("Enter sphere radius to use for surface energy calculation (e.g. 3)\n");
scanf("%d",&rade);
printf("%d \n",rade);
rflag=0; /* always initialize system */
nsph=maketemp(rade);
printf("nsph is %d \n",nsph);
if(rflag==0){
sysinit(ph1id,ph2id);
}
else{
sysscan(ph1id,ph2id);
}
i=0;
rhnow=rhcalc(ph1id);
while((rhnow<rhtarget)&&(i<MAXCYC)){
printf("Now: %f Target: %f \n",rhnow,rhtarget);
i+=1;
printf("Cycle: %d \n",i);
keepgo=movepix(natonce,ph1id,ph2id);
/* If equilibrium is reached, then return to calling routine */
if(keepgo==1){
return;
}
curvsum1=0;
curvsum2=0;
pixsum1=0;
pixsum2=0;
/* Determine average curvatures for phases 1 and 2 */
for(j=0;j<=nsph;j++){
pixsum1+=nsolid[j];
curvsum1+=(j*nsolid[j]);
pixsum2+=nair[j];
curvsum2+=(j*nair[j]);
}
avecurv1=(float)curvsum1/(float)pixsum1;
avecurv2=(float)curvsum2/(float)pixsum2;
printf("Ave. solid curvature: %f \n",avecurv1);
printf("Ave. air curvature: %f \n",avecurv2);
rhnow=rhcalc(ph1id);
}
}
/* routine to read in microstructure from a file */
/* Calls no other routines */
/* Called by main routine */
void readmic()
{
FILE *infile;
char filen[80];
int iout,i1,i2,i3;
long int porcnt,totcnt;
porcnt=totcnt=0;
printf("Enter name of file to read in \n");
scanf("%s",filen);
infile=fopen(filen,"r");
for(i3=0;i3<SYSSIZE;i3++){
for(i2=0;i2<SYSSIZE;i2++){
for(i1=0;i1<SYSSIZE;i1++){
fscanf(infile,"%d",&iout);
totcnt+=1;
phase [i1] [i2] [i3]=iout;
if(iout==0){
porcnt+=1;
}
}
}
}
printf("Count for porosity is %ld out of %ld \n",porcnt,totcnt);
fclose(infile);
}
/* routine to output microstructure to file */
/* Calls no other routines */
/* Called by main routine */
void savemic()
{
FILE *outfile;
int iout,i1,i2,i3;
long int porcnt,totcnt;
char fileout[80];
porcnt=totcnt=0;
printf("Enter name of file to write to \n");
scanf("%s",fileout);
outfile=fopen(fileout,"w");
for(i3=0;i3<SYSSIZE;i3++){
for(i2=0;i2<SYSSIZE;i2++){
for(i1=0;i1<SYSSIZE;i1++){
totcnt+=1;
iout=phase [i1] [i2] [i3];
if(iout==0){
porcnt+=1;
}
fprintf(outfile,"%d\n",iout);
}
}
}
printf("Count for porosity is %ld out of %ld \n",porcnt,totcnt);
fclose(outfile);
}
/* Calls phcount, runsint, readmic, and savemic */
main()
{
int nseed,menuch;
printf("Enter random number seed (integer < 0): \n");
scanf("%d",&nseed);
printf("%d \n",nseed);
seed=(&nseed);
menuch=6;
/* Present menu and obtain user choice */
while(menuch!=1){
printf("Enter choice: \n");
printf("1) Exit program \n");
printf("2) Read in microstructure from file \n");
printf("3) Measure phase fractions \n");
printf("4) Perform sintering algorithm \n");
printf("5) Output resultant microstructure to file \n");
scanf("%d",&menuch);
printf("%d \n",menuch);
switch(menuch){
case 2:
readmic();
break;
case 3:
phcount();
break;
case 4:
runsint();
break;
case 5:
savemic();
break;
default:
break;
}
}
}