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

Validation of Mean Curvature

Browse Source
This commit is contained in:
James E McClure
2013-12-03 16:36:45 -05:00
parent 67fd795b1a
commit 9dd97d7f0f
3 changed files with 344 additions and 37 deletions
Download Patch File Download Diff File Expand all files Collapse all files

44
include/pmmc.h
View File

@@ -3221,13 +3221,13 @@ inline void ComputeAreasPMMC(IntArray &cubeList, int start, int finish,
}
}
//--------------------------------------------------------------------------------------------------------
inline void pmmc_ConstructLocalCube(DoubleArray &SignDist, DoubleArray &Phase, double fluid_isovalue, double solid_isovalue,
inline void pmmc_ConstructLocalCube(DoubleArray &SignDist, DoubleArray &Phase, double solid_isovalue, double fluid_isovalue,
DTMutableList<Point> &nw_pts, IntArray &nw_tris, DoubleArray &values,
DTMutableList<Point> &ns_pts, IntArray &ns_tris,
DTMutableList<Point> &ws_pts, IntArray &ws_tris,
DTMutableList<Point> &local_nws_pts, DTMutableList<Point> &nws_pts, IntArray &nws_seg,
DTMutableList<Point> &local_sol_pts, IntArray &local_sol_tris,
int &n_local_sol_tris, int &n_local_sol_pts, int &n_nw_pts, int n_nw_tris,
int &n_local_sol_tris, int &n_local_sol_pts, int &n_nw_pts, int &n_nw_tris,
int &n_ws_pts, int &n_ws_tris, int &n_ns_tris, int &n_ns_pts,
int &n_local_nws_pts, int &n_nws_pts, int &n_nws_seg,
int i, int j, int k, int Nx, int Ny, int Nz)
@@ -3253,14 +3253,14 @@ inline void pmmc_ConstructLocalCube(DoubleArray &SignDist, DoubleArray &Phase, d
// find the local solid surface using the regular Marching Cubes algorithm
SolidMarchingCubes(SignDist,0.0,Phase,fluid_isovalue,i,j,k,Nx,Ny,Nz,local_sol_pts,n_local_sol_pts,
local_sol_tris,n_local_sol_tris,values);
local_sol_tris,n_local_sol_tris,values);
/////////////////////////////////////////
//////// TRIM THE SOLID SURFACE /////////
/////////////////////////////////////////
TRIM(local_sol_pts, n_local_sol_pts, fluid_isovalue,local_sol_tris, n_local_sol_tris,
ns_pts, n_ns_pts, ns_tris, n_ns_tris, ws_pts, n_ws_pts,
ws_tris, n_ws_tris, values, local_nws_pts, n_local_nws_pts);
ns_pts, n_ns_pts, ns_tris, n_ns_tris, ws_pts, n_ws_pts,
ws_tris, n_ws_tris, values, local_nws_pts, n_local_nws_pts);
/////////////////////////////////////////
//////// WRITE COMMON LINE POINTS ///////
@@ -3274,8 +3274,8 @@ inline void pmmc_ConstructLocalCube(DoubleArray &SignDist, DoubleArray &Phase, d
for (p=0; p<n_local_nws_pts-1; p++){
P = local_nws_pts(p);
if ( P.x == 1.0*i || P.x ==1.0*(i+1)||
P.y == 1.0*j || P.y == 1.0*(j+1) ||
P.z == 1.0*k || P.z == 1.0*(k+1) ){
P.y == 1.0*j || P.y == 1.0*(j+1) ||
P.z == 1.0*k || P.z == 1.0*(k+1) ){
if (p%2 == 0){
// even points
// Swap the pair of points
@@ -3295,7 +3295,7 @@ inline void pmmc_ConstructLocalCube(DoubleArray &SignDist, DoubleArray &Phase, d
}
// guarantee exit from the loop
}
}
}
// Two common curve points per triangle
// 0-(1=2)-(3=4)-...
for (p=1; p<n_local_nws_pts-1; p+=2){
@@ -3378,7 +3378,7 @@ inline void pmmc_MeshGradient(DoubleArray &f, DoubleArray &fx, DoubleArray &fy,
}
}
//--------------------------------------------------------------------------------------------------------
inline void pmmc_MeshCurvature(DoubleArray &f, DoubleArray &MeanCurvature, DoubleArray &GaussCurvature,
inline void pmmc_MeshCurvature(DoubleArray &f, DoubleArray &MeanCurvature, DoubleArray &GaussCurvature,
int Nx, int Ny, int Nz)
{
// Mesh spacing is taken to be one to simplify the calculation
@@ -3397,15 +3397,17 @@ inline void pmmc_MeshCurvature(DoubleArray &f, DoubleArray &MeanCurvature, Doub
fxx = f(i+1,j,k) - 2.0*f(i,j,k) + f(i-1,j,k);
fyy = f(i,j+1,k) - 2.0*f(i,j,k) + f(i,j-1,k);
fzz = f(i,j,k+1) - 2.0*f(i,j,k) + f(i,j,k-1);
fxy = 0.25*(f(i+1,j+1,k) - f(i-1,j-1,k) - f(i-1,j+1,k) + f(i+1,j-1,k));
fxz = 0.25*(f(i+1,j,k+1) - f(i-1,j,k-1) - f(i-1,j,k+1) + f(i+1,j,k-1));
fyz = 0.25*(f(i,j+1,k+1) - f(i,j-1,k-1) - f(i,j-1,k+1) + f(i,j+1,k-1));
fxy = 0.25*(f(i+1,j+1,k) - f(i+1,j-1,k) - f(i-1,j+1,k) + f(i-1,j-1,k));
fxz = 0.25*(f(i+1,j,k+1) - f(i+1,j,k-1) - f(i-1,j,k+1) + f(i-1,j,k-1));
fyz = 0.25*(f(i,j+1,k+1) - f(i,j+1,k-1) - f(i,j-1,k+1) + f(i,j-1,k-1));
// Evaluate the Mean Curvature
denominator = sqrt(pow(fx*fx + fy*fy + fz*fz,3));
denominator = pow(sqrt(fx*fx + fy*fy + fz*fz),3);
if (denominator == 0.0) denominator = 1.0;
MeanCurvature(i,j,k)=(1.0/denominator)*((fyy+fzz)*fx*fx + (fxx+fzz)*fy*fy + (fxx+fyy)*fz*fz
-2.0*fx*fy*fxy - 2.0*fx*fz*fxz - 2.0*fy*fz*fyz);
// Evaluate the Gaussian Curvature
denominator = pow(fx*fx + fy*fy + fz*fz,2);
if (denominator == 0.0) denominator = 1.0;
GaussCurvature(i,j,k) = (1.0/denominator)*(fx*fx*(fyy*fzz-fyz*fyz) + fy*fy*(fxx*fzz-fxz*fxz) + fz*fz*(fxx*fyy-fxy*fxy)
+2.0*(fx*fy*(fxz*fyz-fxy*fzz) + fy*fz*(fxy*fxz-fyz*fxx)
+ fx*fz*(fxy*fyz-fxz*fyy)));
@@ -3460,10 +3462,12 @@ inline double pmmc_CubeSurfaceArea(DTMutableList<Point> &Points, IntArray &Trian
}
//--------------------------------------------------------------------------------------------------------
inline double pmmc_CubeSurfaceInterpValue(DoubleArray &CubeValues, DTMutableList<Point> &Points, IntArray &Triangles,
DoubleArray &SurfaceValues, int npts, int ntris)
DoubleArray &SurfaceValues, int i, int j, int k, int npts, int ntris)
{
Point A,B,C;
int p;
double vA,vB,vC;
double x,y,z;
double s,s1,s2,s3,temp;
double a,b,c,d,e,f,g,h;
double integral;
@@ -3479,9 +3483,12 @@ inline double pmmc_CubeSurfaceInterpValue(DoubleArray &CubeValues, DTMutableList
g = CubeValues(0,1,1)-a-c-d;
h = CubeValues(1,1,1)-a-b-c-d-e-f-g;
for (int i=0; i<npts; i++){
A = Points(i);
SurfaceValues(i) = a + b*A.x + c*A.y+d*A.z + e*A.x*A.y + f*A.x*A.z + g*A.y*A.z + h*A.x*A.y*A.z;
for (p=0; p<npts; p++){
A = Points(p);
x = A.x-1.0*i;
y = A.y-1.0*j;
z = A.z-1.0*k;
SurfaceValues(p) = a + b*x + c*y+d*z + e*x*y + f*x*z + g*y*z + h*x*y*z;
}
integral = 0.0;
@@ -3509,6 +3516,7 @@ inline double pmmc_CubeCurveInterpValue(DoubleArray &CubeValues, DoubleArray &Cu
int p;
Point A,B;
double vA,vB;
double x,y,z;
double s,s1,s2,s3,temp;
double a,b,c,d,e,f,g,h;
double integral;
@@ -3525,7 +3533,7 @@ inline double pmmc_CubeCurveInterpValue(DoubleArray &CubeValues, DoubleArray &Cu
g = CubeValues(0,1,1)-a-c-d;
h = CubeValues(1,1,1)-a-b-c-d-e-f-g;
for (int p=0; p<npts; p++){
for (p=0; p<npts; p++){
A = Points(p);
CurveValues(p) = a + b*A.x + c*A.y+d*A.z + e*A.x*A.y + f*A.x*A.z + g*A.y*A.z + h*A.x*A.y*A.z;
}

283
tests/TestCylinderAreas.cpp Normal file
View File

@@ -0,0 +1,283 @@
#include <iostream>
#include <math.h>
#include "pmmc.h"
//#include "PointList.h"
//#include "Array.h"
#define RADIUS 15
#define HEIGHT 15.5
#define N 60
#define PI 3.14159
int main (int argc, char *argv[])
{
// printf("Radius = %s \n,"RADIUS);
int SIZE = N*N*N;
int Nx,Ny,Nz;
Nx = Ny = Nz = N;
int i,j,k,p,q,r;
// double *Solid; // cylinder
// double *Phase; // region of the cylinder
// Solid = new double [SIZE];
// Phase = new double [SIZE];
DoubleArray SignDist(Nx,Ny,Nz);
DoubleArray Phase(Nx,Ny,Nz);
double fluid_isovalue = 0.0;
double solid_isovalue = 0.0;
/* ****************************************************************
VARIABLES FOR THE PMMC ALGORITHM
****************************************************************** */
//...........................................................................
// Averaging variables
//...........................................................................
double awn,ans,aws,lwns,nwp_volume;
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 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);
// IntArray store;
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;
// 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;
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;
Cx = Cy = Cz = N*0.51;
for (k=0; k<N; k++){
for (j=0; j<N; j++){
for (i=0; i<N; i++){
dist1 = sqrt((i-Cx)*(i-Cx)+(j-Cy)*(j-Cy)) - RADIUS;
dist2 = fabs(Cz-k)-HEIGHT;
// printf("distances = %f, %f \n",dist1,dist2);
//Solid.data[k*Nx*Ny+j*Nx+i] = dist1;
//Phase[k*Nx*Ny+j*Nx+i] = dist2;
SignDist(i,j,k) = -dist1;
Phase(i,j,k) = dist2;
}
}
}
FILE *STRIS;
STRIS = fopen("solid-triangles.out","w");
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");
// End of the loop to set the values
awn = aws = ans = lwns = 0.0;
nwp_volume = 0.0;
As = 0.0;
for (c=0;c<ncubes;c++){
// Get cube from the list
i = cubeList(0,c);
j = cubeList(1,c);
k = cubeList(2,c);
for (p=0;p<8;p++){
if ( Phase(i+cube[p][0],j+cube[p][1],k+cube[p][2]) > 0
&& SignDist(i+cube[p][0],j+cube[p][1],k+cube[p][2]) > 0 ){
nwp_volume += 0.125;
}
}
// 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 Interfacial Areas, Common Line length for blob p
// nw surface
double temp;
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));
// Compute length of sides (assume dx=dy=dz)
s1 = sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)+(A.z-B.z)*(A.z-B.z));
s2 = sqrt((A.x-C.x)*(A.x-C.x)+(A.y-C.y)*(A.y-C.y)+(A.z-C.z)*(A.z-C.z));
s3 = sqrt((B.x-C.x)*(B.x-C.x)+(B.y-C.y)*(B.y-C.y)+(B.z-C.z)*(B.z-C.z));
s = 0.5*(s1+s2+s3);
temp = s*(s-s1)*(s-s2)*(s-s3);
if (temp > 0.0) awn += sqrt(temp);
}
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));
// Compute length of sides (assume dx=dy=dz)
s1 = sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)+(A.z-B.z)*(A.z-B.z));
s2 = sqrt((A.x-C.x)*(A.x-C.x)+(A.y-C.y)*(A.y-C.y)+(A.z-C.z)*(A.z-C.z));
s3 = sqrt((B.x-C.x)*(B.x-C.x)+(B.y-C.y)*(B.y-C.y)+(B.z-C.z)*(B.z-C.z));
s = 0.5*(s1+s2+s3);
//ans=ans+sqrt(s*(s-s1)*(s-s2)*(s-s3));
temp = s*(s-s1)*(s-s2)*(s-s3);
if (temp > 0.0) ans += sqrt(temp);
}
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));
// Compute length of sides (assume dx=dy=dz)
s1 = sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)+(A.z-B.z)*(A.z-B.z));
s2 = sqrt((A.x-C.x)*(A.x-C.x)+(A.y-C.y)*(A.y-C.y)+(A.z-C.z)*(A.z-C.z));
s3 = sqrt((B.x-C.x)*(B.x-C.x)+(B.y-C.y)*(B.y-C.y)+(B.z-C.z)*(B.z-C.z));
s = 0.5*(s1+s2+s3);
//aws=aws+sqrt(s*(s-s1)*(s-s2)*(s-s3));
temp = s*(s-s1)*(s-s2)*(s-s3);
if (temp > 0.0) aws += sqrt(temp);
}
for (r=0;r<n_local_sol_tris;r++){
A = local_sol_pts(local_sol_tris(0,r));
B = local_sol_pts(local_sol_tris(1,r));
C = local_sol_pts(local_sol_tris(2,r));
// Compute length of sides (assume dx=dy=dz)
s1 = sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)+(A.z-B.z)*(A.z-B.z));
s2 = sqrt((A.x-C.x)*(A.x-C.x)+(A.y-C.y)*(A.y-C.y)+(A.z-C.z)*(A.z-C.z));
s3 = sqrt((B.x-C.x)*(B.x-C.x)+(B.y-C.y)*(B.y-C.y)+(B.z-C.z)*(B.z-C.z));
s = 0.5*(s1+s2+s3);
//aws=aws+sqrt(s*(s-s1)*(s-s2)*(s-s3));
temp = s*(s-s1)*(s-s2)*(s-s3);
if (temp > 0.0) As += sqrt(temp);
}
for (p=0; p < n_local_nws_pts-1; p++){
// Extract the line segment
A = local_nws_pts(p);
B = local_nws_pts(p+1);
// Compute the length of the segment
s = sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)+(A.z-B.z)*(A.z-B.z));
// Add the length to the common line
lwns += s;
}
//.......................................................................................
// 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);
}
//.......................................................................................
for (r=0;r<n_local_sol_tris;r++){
A = local_sol_pts(local_sol_tris(0,r));
B = local_sol_pts(local_sol_tris(1,r));
C = local_sol_pts(local_sol_tris(2,r));
fprintf(STRIS,"%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);
}
//*******************************************************************
// Reset the triangle counts to zero
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;
// n_ns_tris_beg = 0;//n_ns_tris;
// n_ws_tris_beg = 0;//n_ws_tris;
// n_nws_seg_beg = n_nws_seg;
//*******************************************************************
}
fclose(WN_TRIS);
fclose(NS_TRIS);
fclose(WS_TRIS);
fclose(WNS_PTS);
fclose(STRIS);
printf("-------------------------------- \n");
printf("NWP volume = %f \n", nwp_volume);
printf("Area wn = %f, Analytical = %f \n", awn,2*PI*RADIUS*RADIUS);
printf("Area ns = %f, Analytical = %f \n", ans, 2*PI*RADIUS*(N-2)-4*PI*RADIUS*HEIGHT);
printf("Area ws = %f, Analytical = %f \n", aws, 4*PI*RADIUS*HEIGHT);
printf("Area s = %f, Analytical = %f \n", As, 2*PI*RADIUS*(N-2));
printf("Length wns = %f, Analytical = %f \n", lwns, 4*PI*RADIUS);
printf("-------------------------------- \n");
//.........................................................................
/* FILE *PHASE;
PHASE = fopen("Phase.in","wb");
fwrite(Phase,8,SIZE,PHASE);
fclose(PHASE);
FILE *SOLID;
SOLID = fopen("Distance.in","wb");
fwrite(Solid,8,SIZE,SOLID);
fclose(SOLID);
*/
}

54
tests/TestSphereCurvature.cpp
View File

@@ -80,13 +80,24 @@ int main(int argc, char **argv)
DoubleArray CubeValues(2,2,2);
// Compute the signed distance function
SignedDistance(SignDist.data,nspheres,cx,cy,cz,rad,Lx,Ly,Lz,Nx,Ny,Nz,0,0,0,1,1,1);
SignedDistance(Phase.data,nspheres,cx,cy,cz,rad,Lx,Ly,Lz,Nx,Ny,Nz,0,0,0,1,1,1);
pmmc_MeshCurvature(SignDist, MeanCurvature, GaussCurvature, Nx, Ny, Nz);
for (k=0; k<Nz; k++){
for (j=0; j<Ny; j++){
for (i=0; i<Nx; i++){
SignDist(i,j,k) = 100.0;
Phase(i,j,k) = sqrt((1.0*i-0.5*Nx)*(1.0*i-0.5*Nx)+(1.0*j-0.5*Ny)*(1.0*j-0.5*Ny)+(1.0*k-0.5*Nz)*(1.0*k-0.5*Nz))-0.3*Nx;
}
}
}
SignedDistance(SignDist.data,0,cx,cy,cz,rad,Lx,Ly,Lz,Nx,Ny,Nz,0,0,0,1,1,1);
pmmc_MeshCurvature(Phase, MeanCurvature, GaussCurvature, Nx, Ny, Nz);
double wn_curvature_sum = 0.0;
double wn_area_sum = 0.0;
for (int c=0;c<ncubes;c++){
// Get cube from the list
@@ -94,6 +105,22 @@ int main(int argc, char **argv)
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);
// Copy the curvature values for the cube
CubeValues(0,0,0) = MeanCurvature(i,j,k);
CubeValues(1,0,0) = MeanCurvature(i+1,j,k);
@@ -104,26 +131,15 @@ int main(int argc, char **argv)
CubeValues(0,1,1) = MeanCurvature(i,j+1,k+1);
CubeValues(1,1,1) = MeanCurvature(i+1,j+1,k+1);
// Construct the interfaces and common curve
pmmc_ConstructLocalCube(SignDist, Phase, fluid_isovalue, solid_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);
// Interpolate the curvature onto the surface
wn_curvature_sum += pmmc_CubeSurfaceInterpValue(CubeValues, local_sol_pts, local_sol_tris,
wn_curvature, n_local_sol_pts, n_local_sol_tris);
wn_curvature_sum += pmmc_CubeSurfaceInterpValue(CubeValues, nw_pts, nw_tris,
wn_curvature, i, j, k, n_nw_pts, n_nw_tris);
wn_area_sum += pmmc_CubeSurfaceArea(nw_pts, nw_tris, n_nw_tris);
wn_area_sum += pmmc_CubeSurfaceArea(local_sol_pts, local_sol_tris, n_local_sol_tris);
}
printf("Area value = %f \n", wn_area_sum);
printf("Curvature sum = %f \n", wn_curvature_sum);
printf("Curvature value = %f, Analytical = %f \n", wn_curvature_sum/wn_area_sum, 2.0/rad[0]/101 );
printf("Mean Curvature Average = %f, Analytical = %f \n", wn_curvature_sum/wn_area_sum, 2.0/rad[0]/101 );
FILE *CURVATURE;
CURVATURE = fopen("Curvature.dat","wb");
@@ -132,7 +148,7 @@ int main(int argc, char **argv)
FILE *DISTANCE;
DISTANCE = fopen("SignDist.dat","wb");
fwrite(SignDist.data,8,N,DISTANCE);
fwrite(Phase.data,8,N,DISTANCE);
fclose(DISTANCE);
}