Added analysis for old LBM codes
This commit is contained in:
parent
42338802c5
commit
da90f34a76
495
pmmc/AnalyzeLBM.cpp
Normal file
495
pmmc/AnalyzeLBM.cpp
Normal file
|
@ -0,0 +1,495 @@
|
|||
#include <iostream>
|
||||
#include <math.h>
|
||||
#include "pmmc.h"
|
||||
#include "Domain.h"
|
||||
//#include "PointList.h"
|
||||
//#include "Array.h"
|
||||
|
||||
using namespace std;
|
||||
|
||||
int main(int argc, char **argv)
|
||||
{
|
||||
//.......................................................................
|
||||
// printf("Radius = %s \n,"RADIUS);
|
||||
int Nx,Ny,Nz,N;
|
||||
int i,j,k,p,q,r,n;
|
||||
int ReadSize; // Number of entries in the input files that are read
|
||||
int nspheres;
|
||||
double Lx,Ly,Lz;
|
||||
//.......................................................................
|
||||
|
||||
//.......................................................................
|
||||
string PhiName,PressName,VelName;
|
||||
ifstream readnames("Filenames.in");
|
||||
readnames >> PhiName;
|
||||
readnames >> PressName;
|
||||
readnames >> VelName;
|
||||
readnames >> ReadSize;
|
||||
//.......................................................................
|
||||
// Reading the domain information file
|
||||
//.......................................................................
|
||||
ifstream domain("Domain.in");
|
||||
domain >> Nx;
|
||||
domain >> Ny;
|
||||
domain >> Nz;
|
||||
domain >> nspheres;
|
||||
domain >> Lx;
|
||||
domain >> Ly;
|
||||
domain >> Lz;
|
||||
//.......................................................................
|
||||
N = Nx*Ny*Nz;
|
||||
printf("Domain size = %i x %i x %i \n",Nx,Ny,Nz);
|
||||
printf("Domain length = %f x %f x %f \n",Lx,Ly,Lz);
|
||||
|
||||
cout << "PhiName = " << PressName << endl;
|
||||
cout << "ReadSize = " << ReadSize << endl;
|
||||
|
||||
//.......................................................................
|
||||
DoubleArray SignDist(Nx,Ny,Nz);
|
||||
DoubleArray SignDistCopy(Nx+2,Ny+2,Nz+2);
|
||||
DoubleArray Phase(Nx,Ny,Nz);
|
||||
DoubleArray Phase_x(Nx,Ny,Nz);
|
||||
DoubleArray Phase_y(Nx,Ny,Nz);
|
||||
DoubleArray Phase_z(Nx,Ny,Nz);
|
||||
DoubleArray Sx(Nx,Ny,Nz);
|
||||
DoubleArray Sy(Nx,Ny,Nz);
|
||||
DoubleArray Sz(Nx,Ny,Nz);
|
||||
DoubleArray Vel_x(Nx,Ny,Nz);
|
||||
DoubleArray Vel_y(Nx,Ny,Nz);
|
||||
DoubleArray Vel_z(Nx,Ny,Nz);
|
||||
DoubleArray Press(Nx,Ny,Nz);
|
||||
DoubleArray GaussCurvature(Nx,Ny,Nz);
|
||||
DoubleArray MeanCurvature(Nx,Ny,Nz);
|
||||
//.......................................................................
|
||||
|
||||
// Read in sphere pack
|
||||
printf("nspheres =%i \n",nspheres);
|
||||
//.......................................................................
|
||||
double *cx,*cy,*cz,*rad;
|
||||
cx = new double[nspheres];
|
||||
cy = new double[nspheres];
|
||||
cz = new double[nspheres];
|
||||
rad = new double[nspheres];
|
||||
//.......................................................................
|
||||
int *ReadInfo;
|
||||
double *ReadPhi, *ReadPress, *ReadVel;
|
||||
ReadInfo = new int[3*ReadSize];
|
||||
ReadPhi = new double[ReadSize];
|
||||
ReadPress = new double[ReadSize];
|
||||
ReadVel = new double[3*ReadSize];
|
||||
//.......................................................................
|
||||
printf("Reading the sphere packing \n");
|
||||
ReadSpherePacking(nspheres,cx,cy,cz,rad);
|
||||
|
||||
SignedDistance(SignDistCopy.data,nspheres,cx,cy,cz,rad,Lx,Ly,Lz,Nx+2,Ny+2,Nz+2,
|
||||
0,0,0,1,1,1);
|
||||
|
||||
for (k=0; k<Nz; k++){
|
||||
for (j=0; j<Ny; j++){
|
||||
for (i=0; i<Nx; i++){
|
||||
SignDist(i,j,k) = SignDistCopy(i+1,j+1,k+1)-1.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
//.......................................................................
|
||||
// Read all of the output files from Lattice Boltzmann Simulation
|
||||
//.......................................................................
|
||||
printf("Reading data from the Lattice Boltzmann output files \n");
|
||||
int infovalue;
|
||||
double value;
|
||||
ifstream INFO("Info",ios::binary);
|
||||
for (n=0; n< ReadSize; n++){
|
||||
// Read 3D indices from the info file
|
||||
INFO.read((char*) &infovalue, sizeof(infovalue));
|
||||
ReadInfo[3*n] = infovalue;
|
||||
INFO.read((char*) &infovalue, sizeof(infovalue));
|
||||
ReadInfo[3*n+1] = infovalue;
|
||||
INFO.read((char*) &infovalue, sizeof(infovalue));
|
||||
ReadInfo[3*n+2] = infovalue-4; // account for pressure boundary conditions
|
||||
}
|
||||
INFO.close();
|
||||
ifstream PHI(PhiName.c_str(),ios::binary);
|
||||
for (n=0; n< ReadSize; n++){
|
||||
PHI.read((char*) &value, sizeof(value));
|
||||
ReadPhi[n] = value;
|
||||
}
|
||||
PHI.close();
|
||||
ifstream PRESSURE(PressName.c_str(),ios::binary);
|
||||
for (n=0; n< ReadSize; n++){
|
||||
PRESSURE.read((char*) &value, sizeof(value));
|
||||
ReadPress[n] = value;
|
||||
}
|
||||
PRESSURE.close();
|
||||
|
||||
#ifdef ReadVelocityFromFile
|
||||
ifstream VELOCITY(VelName.c_str(),ios::binary);
|
||||
for (n=0; n< ReadSize; n++){
|
||||
VELOCITY.read((char*) &value, sizeof(value));
|
||||
ReadVel[3*n] = value;
|
||||
VELOCITY.read((char*) &value, sizeof(value));
|
||||
ReadVel[3*n+1] = value;
|
||||
VELOCITY.read((char*) &value, sizeof(value));
|
||||
ReadVel[3*n+2] = value;
|
||||
}
|
||||
VELOCITY.close();
|
||||
#else
|
||||
value = 0.f;
|
||||
for (n=0; n< ReadSize; n++){
|
||||
ReadVel[3*n] = value;
|
||||
ReadVel[3*n+1] = value;
|
||||
ReadVel[3*n+2] = value;
|
||||
}
|
||||
#endif
|
||||
//.......................................................................
|
||||
|
||||
//.......................................................................
|
||||
// Reconstruct the global meshes using the info file
|
||||
//.......................................................................
|
||||
for (int idx=0; idx< ReadSize; idx++){
|
||||
// 3D indices
|
||||
i = ReadInfo[3*idx];
|
||||
j = ReadInfo[3*idx+1];
|
||||
k = ReadInfo[3*idx+2];
|
||||
// printf("%i, %i, %i \n",i,j,k);
|
||||
// Update the values
|
||||
if (k>0){
|
||||
Phase(i,j,k) = ReadPhi[idx];
|
||||
Press(i,j,k) = ReadPress[idx];
|
||||
Vel_x(i,j,k) = ReadVel[3*idx];
|
||||
Vel_y(i,j,k) = ReadVel[3*idx+1];
|
||||
Vel_z(i,j,k) = ReadVel[3*idx+2];
|
||||
}
|
||||
}
|
||||
//.......................................................................
|
||||
|
||||
//.......................................................................
|
||||
double fluid_isovalue = 0.0;
|
||||
double solid_isovalue = 0.0;
|
||||
//.......................................................................
|
||||
|
||||
/* //.......................................................................
|
||||
double *cx,*cy,*cz,*rad;
|
||||
cx = new double[nspheres];
|
||||
cy = new double[nspheres];
|
||||
cz = new double[nspheres];
|
||||
rad = new double[nspheres];
|
||||
//...............................
|
||||
printf("Reading the sphere packing \n");
|
||||
ReadSpherePacking(nspheres,cx,cy,cz,rad);
|
||||
//.......................................................................
|
||||
//.......................................................................
|
||||
// Compute the signed distance function for the sphere packing
|
||||
SignedDistance(SignDist.data,nspheres,cx,cy,cz,rad,Lx,Ly,Lz,Nx,Ny,Nz,0,0,0,1,1,1);
|
||||
*/ //.......................................................................
|
||||
|
||||
/* ****************************************************************
|
||||
VARIABLES FOR THE PMMC ALGORITHM
|
||||
****************************************************************** */
|
||||
//...........................................................................
|
||||
// Averaging variables
|
||||
//...........................................................................
|
||||
double awn,ans,aws,lwns,nwp_volume;
|
||||
double sw,vol_n,vol_w,paw,pan;
|
||||
double efawns,Jwn;
|
||||
double As;
|
||||
double dEs,dAwn,dAns; // Global surface energy (calculated by rank=0)
|
||||
double awn_global,ans_global,aws_global,lwns_global,nwp_volume_global;
|
||||
double As_global;
|
||||
// bool add=1; // Set to false if any corners contain nw-phase ( F > fluid_isovalue)
|
||||
|
||||
int n_nw_pts=0,n_ns_pts=0,n_ws_pts=0,n_nws_pts=0, map=0;
|
||||
int n_nw_tris=0, n_ns_tris=0, n_ws_tris=0, n_nws_seg=0;
|
||||
|
||||
double s,s1,s2,s3; // Triangle sides (lengths)
|
||||
Point A,B,C,P;
|
||||
// double area;
|
||||
int cube[8][3] = {{0,0,0},{1,0,0},{0,1,0},{1,1,0},{0,0,1},{1,0,1},{0,1,1},{1,1,1}}; // cube corners
|
||||
// int count_in=0,count_out=0;
|
||||
// int nodx,nody,nodz;
|
||||
// initialize lists for vertices for surfaces, common line
|
||||
DTMutableList<Point> nw_pts(20);
|
||||
DTMutableList<Point> ns_pts(20);
|
||||
DTMutableList<Point> ws_pts(20);
|
||||
DTMutableList<Point> nws_pts(20);
|
||||
// initialize triangle lists for surfaces
|
||||
IntArray nw_tris(3,20);
|
||||
IntArray ns_tris(3,20);
|
||||
IntArray ws_tris(3,20);
|
||||
// initialize list for line segments
|
||||
IntArray nws_seg(2,20);
|
||||
DTMutableList<Point> tmp(20);
|
||||
|
||||
// Initialize arrays for local solid surface
|
||||
DTMutableList<Point> local_sol_pts(20);
|
||||
int n_local_sol_pts = 0;
|
||||
IntArray local_sol_tris(3,18);
|
||||
int n_local_sol_tris;
|
||||
DoubleArray values(20);
|
||||
DTMutableList<Point> local_nws_pts(20);
|
||||
int n_local_nws_pts;
|
||||
|
||||
DoubleArray CubeValues(2,2,2);
|
||||
DoubleArray ContactAngle(20);
|
||||
DoubleArray Curvature(20);
|
||||
DoubleArray InterfaceSpeed(20);
|
||||
DoubleArray NormalVector(60);
|
||||
DoubleArray van(3);
|
||||
DoubleArray vaw(3);
|
||||
DoubleArray vawn(3);
|
||||
DoubleArray Gwn(6);
|
||||
DoubleArray Gns(6);
|
||||
DoubleArray Gws(6);
|
||||
|
||||
double iVol = 1.0/Nx/Ny/Nz;
|
||||
|
||||
int c;
|
||||
//...........................................................................
|
||||
int ncubes = (Nx-2)*(Ny-2)*(Nz-2); // Exclude the "upper" halo
|
||||
IntArray cubeList(3,ncubes);
|
||||
pmmc_CubeListFromMesh(cubeList, ncubes, Nx, Ny, Nz);
|
||||
//...........................................................................
|
||||
double Cx,Cy,Cz;
|
||||
double dist1,dist2;
|
||||
// Extra copies of phase indicator needed to compute time derivatives on CPU
|
||||
DoubleArray Phase_tminus(Nx,Ny,Nz);
|
||||
DoubleArray Phase_tplus(Nx,Ny,Nz);
|
||||
DoubleArray dPdt(Nx,Ny,Nz);
|
||||
for (k=0; k<Nz; k++){
|
||||
for (j=0; j<Ny; j++){
|
||||
for (i=0; i<Nx; i++){
|
||||
Phase_tminus(i,j,k) = Phase(i,j,k);
|
||||
Phase_tplus(i,j,k) = Phase(i,j,k);
|
||||
dPdt(i,j,k) = 0.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
printf("Computing TCAT Averages \n");
|
||||
printf("--------------------------------------------------------------------------------------\n");
|
||||
printf("sw pw pn vw[x, y, z] vn[x, y, z] "); // Volume averages
|
||||
printf("awn ans aws Jwn lwns efawns "); // Interface and common curve averages
|
||||
printf("Gwn [xx, yy, zz, xy, xz, yz] "); // Orientation tensors
|
||||
printf("Gws [xx, yy, zz, xy, xz, yz] ");
|
||||
printf("Gns [xx, yy, zz, xy, xz, yz] \n");
|
||||
printf("--------------------------------------------------------------------------------------\n");
|
||||
|
||||
awn = aws = ans = lwns = 0.0;
|
||||
nwp_volume = 0.0;
|
||||
As = 0.0;
|
||||
Jwn = 0.0;
|
||||
efawns = 0.0;
|
||||
// Compute phase averages
|
||||
pan = paw = 0.0;
|
||||
vaw(0) = vaw(1) = vaw(2) = 0.0;
|
||||
van(0) = van(1) = van(2) = 0.0;
|
||||
vawn(0) = vawn(1) = vawn(2) = 0.0;
|
||||
Gwn(0) = Gwn(1) = Gwn(2) = 0.0;
|
||||
Gwn(3) = Gwn(4) = Gwn(5) = 0.0;
|
||||
Gws(0) = Gws(1) = Gws(2) = 0.0;
|
||||
Gws(3) = Gws(4) = Gws(5) = 0.0;
|
||||
Gns(0) = Gns(1) = Gns(2) = 0.0;
|
||||
Gns(3) = Gns(4) = Gns(5) = 0.0;
|
||||
vol_w = vol_n =0.0;
|
||||
|
||||
// Calculate the time derivative of the phase indicator field
|
||||
for (int n=0; n<Nx*Ny*Nz; n++) dPdt(n) = 0.5*(Phase_tplus(n) - Phase_tminus(n));
|
||||
|
||||
pmmc_MeshGradient(Phase,Phase_x,Phase_y,Phase_z,Nx,Ny,Nz);
|
||||
pmmc_MeshGradient(SignDist,Sx,Sy,Sz,Nx,Ny,Nz);
|
||||
pmmc_MeshCurvature(Phase, MeanCurvature, GaussCurvature, Nx, Ny, Nz);
|
||||
|
||||
/// Compute volume averages
|
||||
int total = 0;
|
||||
// Compute phase averages
|
||||
nwp_volume = 0.0;
|
||||
pan = paw = 0.0;
|
||||
vaw(0) = vaw(1) = vaw(2) = 0.0;
|
||||
vol_w = vol_n =0.0;
|
||||
for (k=0; k<Nz; k++){
|
||||
for (j=0; j<Ny; j++){
|
||||
for (i=0; i<Nx; i++){
|
||||
if ( SignDist(i,j,k) > 0 ){
|
||||
total++;
|
||||
// 1-D index for this cube corner
|
||||
n = i + j*Nx + k*Nx*Ny;
|
||||
|
||||
// Compute the non-wetting phase volume contribution
|
||||
if ( Phase(i,j,k) > 0.0 )
|
||||
nwp_volume += 1.0;
|
||||
|
||||
// volume averages over the non-wetting phase
|
||||
if ( Phase(i,j,k) > 0.9999 ){
|
||||
// volume the excludes the interfacial region
|
||||
vol_n += 1.0;
|
||||
// pressure
|
||||
pan += Press(i,j,k);
|
||||
// velocity
|
||||
van(0) += Vel_x(i,j,k);
|
||||
van(1) += Vel_y(i,j,k);
|
||||
van(2) += Vel_z(i,j,k);
|
||||
}
|
||||
|
||||
// volume averages over the wetting phase
|
||||
if ( Phase(i,j,k) < -0.9999 ){
|
||||
// volume the excludes the interfacial region
|
||||
vol_w += 1.0;
|
||||
// pressure
|
||||
paw += Press(i,j,k);
|
||||
// velocity
|
||||
vaw(0) += Vel_x(i,j,k);
|
||||
vaw(1) += Vel_y(i,j,k);
|
||||
vaw(2) += Vel_z(i,j,k);
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// End of the loop to set the values
|
||||
awn = aws = ans = lwns = 0.0;
|
||||
As = 0.0;
|
||||
Gwn(0) = Gwn(1) = Gwn(2) = Gwn(3) = Gwn(4) = Gwn(5) = 0.0;
|
||||
Gws(0) = Gws(1) = Gws(2) = Gws(3) = Gws(4) = Gws(5) = 0.0;
|
||||
Gns(0) = Gns(1) = Gns(2) = Gns(3) = Gns(4) = Gns(5) = 0.0;
|
||||
|
||||
FILE *WN_TRIS;
|
||||
WN_TRIS = fopen("wn-tris.out","w");
|
||||
|
||||
FILE *NS_TRIS;
|
||||
NS_TRIS = fopen("ns-tris.out","w");
|
||||
|
||||
FILE *WS_TRIS;
|
||||
WS_TRIS = fopen("ws-tris.out","w");
|
||||
|
||||
FILE *WNS_PTS;
|
||||
WNS_PTS = fopen("wns-pts.out","w");
|
||||
|
||||
for (c=0;c<ncubes;c++){
|
||||
// Get cube from the list
|
||||
i = cubeList(0,c);
|
||||
j = cubeList(1,c);
|
||||
k = cubeList(2,c);
|
||||
|
||||
// Run PMMC
|
||||
n_local_sol_tris = 0;
|
||||
n_local_sol_pts = 0;
|
||||
n_local_nws_pts = 0;
|
||||
n_nw_pts=0,n_ns_pts=0,n_ws_pts=0,n_nws_pts=0, map=0;
|
||||
n_nw_tris=0, n_ns_tris=0, n_ws_tris=0, n_nws_seg=0;
|
||||
|
||||
// Construct the interfaces and common curve
|
||||
pmmc_ConstructLocalCube(SignDist, Phase, solid_isovalue, fluid_isovalue,
|
||||
nw_pts, nw_tris, values, ns_pts, ns_tris, ws_pts, ws_tris,
|
||||
local_nws_pts, nws_pts, nws_seg, local_sol_pts, local_sol_tris,
|
||||
n_local_sol_tris, n_local_sol_pts, n_nw_pts, n_nw_tris,
|
||||
n_ws_pts, n_ws_tris, n_ns_tris, n_ns_pts, n_local_nws_pts, n_nws_pts, n_nws_seg,
|
||||
i, j, k, Nx, Ny, Nz);
|
||||
|
||||
// Compute the velocity of the wn interface
|
||||
pmmc_InterfaceSpeed(dPdt, Phase_x, Phase_y, Phase_z, CubeValues, nw_pts, nw_tris,
|
||||
NormalVector, InterfaceSpeed, vawn, i, j, k, n_nw_pts, n_nw_tris);
|
||||
|
||||
// pmmc_InterfaceSpeed(dPdt, Phase_x, Phase_y, Phase_z, CubeValues, nw_pts, nw_tris,
|
||||
// NormalVector, InterfaceSpeed, vawn, i, j, k, n_nw_pts, n_nw_tris);
|
||||
|
||||
|
||||
// Compute the average contact angle
|
||||
efawns += pmmc_CubeContactAngle(CubeValues,ContactAngle,Phase_x,Phase_y,Phase_z,Sx,Sy,Sz,
|
||||
local_nws_pts,i,j,k,n_local_nws_pts);
|
||||
|
||||
// Compute the curvature of the wn interface
|
||||
Jwn += pmmc_CubeSurfaceInterpValue(CubeValues, MeanCurvature, nw_pts, nw_tris,
|
||||
Curvature, i, j, k, n_nw_pts, n_nw_tris);
|
||||
|
||||
// Compute the surface orientation and the interfacial area
|
||||
awn += pmmc_CubeSurfaceOrientation(Gwn,nw_pts,nw_tris,n_nw_tris);
|
||||
ans += pmmc_CubeSurfaceOrientation(Gns,ns_pts,ns_tris,n_ns_tris);
|
||||
aws += pmmc_CubeSurfaceOrientation(Gws,ws_pts,ws_tris,n_ws_tris);
|
||||
|
||||
//*******************************************************************
|
||||
// Compute the Interfacial Areas, Common Line length
|
||||
As += pmmc_CubeSurfaceArea(local_sol_pts,local_sol_tris,n_local_sol_tris);
|
||||
lwns += pmmc_CubeCurveLength(local_nws_pts,n_local_nws_pts);
|
||||
|
||||
//.......................................................................................
|
||||
// Write the triangle lists to text file
|
||||
for (r=0;r<n_nw_tris;r++){
|
||||
A = nw_pts(nw_tris(0,r));
|
||||
B = nw_pts(nw_tris(1,r));
|
||||
C = nw_pts(nw_tris(2,r));
|
||||
fprintf(WN_TRIS,"%f %f %f %f %f %f %f %f %f \n",A.x,A.y,A.z,B.x,B.y,B.z,C.x,C.y,C.z);
|
||||
}
|
||||
for (r=0;r<n_ws_tris;r++){
|
||||
A = ws_pts(ws_tris(0,r));
|
||||
B = ws_pts(ws_tris(1,r));
|
||||
C = ws_pts(ws_tris(2,r));
|
||||
fprintf(WS_TRIS,"%f %f %f %f %f %f %f %f %f \n",A.x,A.y,A.z,B.x,B.y,B.z,C.x,C.y,C.z);
|
||||
}
|
||||
for (r=0;r<n_ns_tris;r++){
|
||||
A = ns_pts(ns_tris(0,r));
|
||||
B = ns_pts(ns_tris(1,r));
|
||||
C = ns_pts(ns_tris(2,r));
|
||||
fprintf(NS_TRIS,"%f %f %f %f %f %f %f %f %f \n",A.x,A.y,A.z,B.x,B.y,B.z,C.x,C.y,C.z);
|
||||
}
|
||||
for (p=0; p < n_nws_pts; p++){
|
||||
P = nws_pts(p);
|
||||
fprintf(WNS_PTS,"%f %f %f \n",P.x, P.y, P.z);
|
||||
}
|
||||
|
||||
}
|
||||
fclose(WN_TRIS);
|
||||
fclose(NS_TRIS);
|
||||
fclose(WS_TRIS);
|
||||
fclose(WNS_PTS);
|
||||
|
||||
printf("paw = %f \n",paw);
|
||||
printf("vol_w = %f \n",vol_w);
|
||||
printf("efawns = %f \n",efawns);
|
||||
printf("lwns = %f \n",lwns);
|
||||
printf("efawns = %f \n",efawns/lwns);
|
||||
|
||||
paw /= vol_w;
|
||||
pan /= vol_n;
|
||||
for (i=0; i<3; i++) vaw(i) /= vol_w;
|
||||
for (i=0; i<3; i++) van(i) /= vol_n;
|
||||
|
||||
Jwn /= awn;
|
||||
efawns /= lwns;
|
||||
for (i=0; i<3; i++) vawn(i) /= awn;
|
||||
for (i=0; i<6; i++) Gwn(i) /= awn;
|
||||
for (i=0; i<6; i++) Gns(i) /= ans;
|
||||
for (i=0; i<6; i++) Gws(i) /= aws;
|
||||
|
||||
awn = awn*iVol;
|
||||
aws = aws*iVol;
|
||||
ans = ans*iVol;
|
||||
lwns = lwns*iVol;
|
||||
sw = 1.0 - nwp_volume/total;
|
||||
|
||||
printf("%.5g %.5g %.5g ",sw,paw,pan); // saturation and pressure
|
||||
printf("%.5g %.5g %.5g ",vaw(0),vaw(1),vaw(2)); // average velocity of w phase
|
||||
printf("%.5g %.5g %.5g ",van(0),van(1),van(2)); // average velocity of n phase
|
||||
printf("%.5g %.5g %.5g ",awn,ans,aws); // interfacial areas
|
||||
printf("%.5g ",Jwn); // curvature of wn interface
|
||||
printf("%.5g ",lwns); // common curve length
|
||||
printf("%.5g ",efawns); // average contact angle
|
||||
printf("%.5g %.5g %.5g %.5g %.5g %.5g ",
|
||||
Gwn(0),Gwn(1),Gwn(2),Gwn(3),Gwn(4),Gwn(5)); // orientation of wn interface
|
||||
printf("%.5g %.5g %.5g %.5g %.5g %.5g ",
|
||||
Gns(0),Gns(1),Gns(2),Gns(3),Gns(4),Gns(5)); // orientation of ns interface
|
||||
printf("%.5g %.5g %.5g %.5g %.5g %.5g \n",
|
||||
Gws(0),Gws(1),Gws(2),Gws(3),Gws(4),Gws(5)); // orientation of ws interface
|
||||
|
||||
FILE *PHASE;
|
||||
PHASE = fopen("Phase.out","wb");
|
||||
fwrite(Phase.data,8,N,PHASE);
|
||||
fclose(PHASE);
|
||||
|
||||
FILE *SOLID;
|
||||
SOLID = fopen("Distance.out","wb");
|
||||
fwrite(SignDist.data,8,N,SOLID);
|
||||
fclose(SOLID);
|
||||
|
||||
}
|
Loading…
Reference in New Issue
Block a user