From 61cbeb4220cc5a7c97ad7312a1c43d8f7ed5c169 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Sat, 16 Apr 2016 19:27:24 -0400 Subject: [PATCH] Updated lbpm_segmented_pp based on Min (2010) to sovle Eikonal equation --- common/Domain.h | 27 ++++--- tests/lbpm_segmented_pp.cpp | 152 +++++++++++++++++++++++++++++++++++- 2 files changed, 167 insertions(+), 12 deletions(-) diff --git a/common/Domain.h b/common/Domain.h index 80ee9572..5072f8d9 100755 --- a/common/Domain.h +++ b/common/Domain.h @@ -137,7 +137,7 @@ static inline void fgetl( char * str, int num, FILE * stream ) -inline void SSO(DoubleArray &Distance, char *ID, Domain &Dm, int timesteps){ +inline double SSO(DoubleArray &Distance, char *ID, Domain &Dm, 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 @@ -152,6 +152,7 @@ inline void SSO(DoubleArray &Distance, char *ID, Domain &Dm, int timesteps){ int in,jn,kn; double Dqx,Dqy,Dqz,Dx,Dy,Dz,W; double nx,ny,nz,Cqx,Cqy,Cqz,sign,norm; + double TotalVariation; const static int D3Q27[26][3]={{1,0,0},{-1,0,0},{0,1,0},{0,-1,0},{0,0,1},{0,0,-1}, {1,1,0},{-1,-1,0},{1,-1,0},{-1,1,0},{1,0,1},{-1,0,-1},{1,0,-1},{-1,0,1}, @@ -176,6 +177,7 @@ inline void SSO(DoubleArray &Distance, char *ID, Domain &Dm, int timesteps){ // Communicate the halo of values fillData.fill(Distance); + TotalVariation=0.0; // Execute the next timestep for (k=1;k 0.0){ + Dx /= W; + Dy /= W; + Dz /= W; + } norm = sqrt(Dx*Dx+Dy*Dy+Dz*Dz); } else{ norm = 0.0; } Distance(i,j,k) += dt*sign*(1.0 - norm); - + TotalVariation += dt*sign*(1.0 - norm); // Disallow any change in phase // if (Distance(i,j,k)*2.0*(ID[n]-1.0) < 0) Distance(i,j,k) = -Distance(i,j,k); } } } - + TotalVariation /= (Dm.Nx-2)*(Dm.Ny-2)*(Dm.Nz-2); count++; } + + return TotalVariation; } diff --git a/tests/lbpm_segmented_pp.cpp b/tests/lbpm_segmented_pp.cpp index 0440f13e..10bc57c1 100644 --- a/tests/lbpm_segmented_pp.cpp +++ b/tests/lbpm_segmented_pp.cpp @@ -26,6 +26,151 @@ inline void MeanFilter(DoubleArray &Mesh){ } } +inline double minmod(double &a, double &b){ + + double value; + + value = a; + if ( if a*b < 0.0) value=0.0; + else if (fabs(a) > fabs(b) value = b; + + return value; +} + + +inline double Eikonal(DoubleArray &Distance, char *ID, Domain &Dm, 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 i,j,k; + double dt=0.25; + double Dx,Dy,Dz; + double Dxp,Dxm,Dyp,Dym,Dzp,Dzm; + double Dxxp,Dxxm,Dyyp,Dyym,Dzzp,Dzzm; + double sign,norm; + double LocalVar,GlobalVar,LocalMax,GlobalMax; + + int xdim,ydim,zdim; + xdim=Dm.Nx-2; + ydim=Dm.Ny-2; + zdim=Dm.Nz-2; + fillHalo fillData(Dm.Comm, Dm.rank_info,xdim,ydim,zdim,1,1,1,0,1); + + // Arrays to store the second derivatives + DoubleArray Dxx(Nx,Ny,Nz); + DoubleArray Dyy(Nx,Ny,Nz); + DoubleArray Dzz(Nx,Ny,Nz); + + int count = 0; + while (count < timesteps){ + + // Communicate the halo of values + fillData.fill(Distance); + + // Compute second order derivatives + for (k=1;k Dxm) Dx = Dxp - Distance(i,j,k) + 0.5*Dxxp; + else Dx = Distance(i,j,k) - Dxm + 0.5*Dxxm; + + if (Dyp > Dym) Dy = Dyp - Distance(i,j,k) + 0.5*Dyyp; + else Dy = Distance(i,j,k) - Dym + 0.5*Dyym; + + if (Dzp > Dzm) Dz = Dzp - Distance(i,j,k) + 0.5*Dzzp; + else Dz = Distance(i,j,k) - Dzm + 0.5*Dzzm; + } + else{ + if (Dxp < Dxm) Dx = Dxp - Distance(i,j,k) + 0.5*Dxxp; + else Dx = Distance(i,j,k) - Dxm + 0.5*Dxxm; + + if (Dyp < Dym) Dy = Dyp - Distance(i,j,k) + 0.5*Dyyp; + else Dy = Distance(i,j,k) - Dym + 0.5*Dyym; + + if (Dzp < Dzm) Dz = Dzp - Distance(i,j,k) + 0.5*Dzzp; + else Dz = Distance(i,j,k) - Dzm + 0.5*Dzzm; + } + + norm=sqrt(Dx*Dx+Dy*Dy+Dz*Dz); + + Distance(i,j,k) += dt*sign*(1.0 - norm); + LocalVar += dt*sign*(1.0 - norm); + + if (fabs(dt*sign*(1.0 - norm)) > LocalMax) + LocalMax = fabs(dt*sign*(1.0 - norm)); + } + } + } + + MPI_Allreduce(&LocalVar,&GlobalVar,1,MPI_DOUBLE,MPI_SUM,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, Global variation=%f \n",count,GlobalVar); + } + + return GlobalVar; +} + + int main(int argc, char **argv) { // Initialize MPI @@ -146,8 +291,13 @@ 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"); - SSO(Averages.SDs,id,Dm,100); + LocalVar = Eikonal(Averages.SDs,id,Dm,100); + + 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); sprintf(LocalRankFilename,"SignDist.%05i",rank); FILE *DIST = fopen(LocalRankFilename,"wb");