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

area to sphere curvature

Browse Source
This commit is contained in:
James E McClure
2018-09-10 17:43:12 -04:00
parent fea083fc82
commit 64df208d45
1 changed files with 17 additions and 6 deletions
Download Patch File Download Diff File Expand all files Collapse all files

23
tests/TestSphereCurvature.cpp
View File

@@ -32,25 +32,36 @@ int main(int argc, char **argv)
}
}
}
double wn_curvature_sum = 0.0;
double wn_area_sum = 0.0;
printf("Construct local isosurface \n");
DECL sphere;
point V1,V2,V3;
unsigned long int e1,e2,e3;
double s,s1,s2,s3;
double area = 0.f;
for (k=0; k<Nz-1; k++){
for (j=0; j<Ny-1; j++){
for (i=0; i<Nx-1; i++){
sphere.LocalIsosurface(Phase,0.f,i,j,k);
for (unsigned long int idx=0; idx<sphere.TriangleCount; idx++){
unsigned long int edge = sphere.Face(idx);
e1 = sphere.Face(idx);
e2 = sphere.halfedge.next(e1);
e3 = sphere.halfedge.next(e2);
V1 = vertex(sphere.halfedge.v1(e1));
V2 = vertex(sphere.halfedge.v1(e2));
V3 = vertex(sphere.halfedge.v1(e3));
s1 = sqrt((V1.x-V2.x)*(V1.x-V2.x)+(V1.y-V2.y)*(V1.y-V2.y)+(V1.z-V2.z)*(V1.z-V2.z));
s2 = sqrt((V1.x-V3.x)*(V1.x-V3.x)+(V1.y-V3.y)*(V1.y-V3.y)+(V1.z-V3.z)*(V1.z-V3.z));
s3 = sqrt((V2.x-V3.x)*(V2.x-V3.x)+(V2.y-V3.y)*(V2.y-V3.y)+(V2.z-V3.z)*(V2.z-V3.z));
s = 0.5*(s1+s2+s3);
area += sqrt(s*(s-s1)*(s-s2)*(s-s3));
}
}
}
}
// printf("Mean Curvature Average = %f, Analytical = %f \n", wn_curvature_sum/wn_area_sum, 2.0/rad[0]/101 );
printf("Surface area = %f (analytical = %f) \n", area,4*3.14159*0.3*0.3*double(Nx*Nx));
// printf("Mean Curvature Average = %f, Analytical = %f \n", wn_curvature_sum/wn_area_sum, 2.0/rad[0]/101 );
int toReturn = 0;
/* if ( fabs(wn_curvature_sum/wn_area_sum -2.0/rad[0]/101)*rad[0]*101.0*0.5 > 0.01 ){