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

still debugging curvature

Browse Source
This commit is contained in:
James E McClure
2018-09-13 13:20:57 -04:00
parent 96d08f6bd8
commit c2b5605e67
2 changed files with 19 additions and 21 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
analysis/decl.cpp
View File

@@ -388,11 +388,11 @@ double DECL::EdgeAngle(int edge)
angle *= 0.5; // half edge value angle *= 0.5; // half edge value
} }
//1.570796326794897 //1.570796326794897
/* if (fabs(angle) > 0.f){ if (fabs(angle) > 0.f){
printf("%f, %f (Edge=%i, twin=%i): U={%f, %f, %f}, V={%f, %f, %f}\n",angle,dotprod,edge,halfedge.twin(edge),U.x,U.y,U.z,V.x,V.y,V.z); printf("%f, %f (Edge=%i, twin=%i): U={%f, %f, %f}, V={%f, %f, %f}\n",angle,dotprod,edge,halfedge.twin(edge),U.x,U.y,U.z,V.x,V.y,V.z);
printf(" P={%f, %f, %f}, Q={%f, %f, %f}, R={%f, %f, %f} \n",P.x,P.y,P.z,Q.x,Q.y,Q.z,R.x,R.y,R.z); printf(" P={%f, %f, %f}, Q={%f, %f, %f}, R={%f, %f, %f} \n",P.x,P.y,P.z,Q.x,Q.y,Q.z,R.x,R.y,R.z);
} }
*/
return angle; return angle;
} }

36
tests/TestSphereCurvature.cpp
View File

@@ -42,6 +42,23 @@ int main(int argc, char **argv)
Nx+=2; Ny+=2; Nz+=2; Nx+=2; Ny+=2; Nz+=2;
DoubleArray Phase(Nx,Ny,Nz); DoubleArray Phase(Nx,Ny,Nz);
Minkowski plane(Dm);
printf("Set distance map \n");
for (k=0; k<Nz; k++){
for (j=0; j<Ny; j++){
for (i=0; i<Nx; i++){
Phase(i,j,k) = k-0.5*double(Nz+1);
}
}
}
printf("Construct local isosurface \n");
plane.ComputeScalar(Phase,0.f);
printf("Surface area = %f (analytical = %f) \n", plane.Ai,double((Nx-2)*(Ny-2)));
printf("Mean curvature = %f (analytical =0) \n", plane.Ji);
printf("Euler characteristic = %f (analytical = 2.0) \n",plane.Xi);
Minkowski sphere(Dm); Minkowski sphere(Dm);
@@ -61,25 +78,6 @@ int main(int argc, char **argv)
printf("Mean curvature = %f (analytical = %f) \n", sphere.Ji,8*3.14159*0.3*double(Nx)); printf("Mean curvature = %f (analytical = %f) \n", sphere.Ji,8*3.14159*0.3*double(Nx));
printf("Euler characteristic = %f (analytical = 2.0) \n",sphere.Xi); printf("Euler characteristic = %f (analytical = 2.0) \n",sphere.Xi);
Minkowski plane(Dm);
printf("Set distance map \n");
for (k=0; k<Nz; k++){
for (j=0; j<Ny; j++){
for (i=0; i<Nx; i++){
Phase(i,j,k) = k-0.5*double(Nz+1);
}
}
}
printf("Construct local isosurface \n");
plane.ComputeScalar(Phase,0.f);
printf("Surface area = %f (analytical = %f) \n", plane.Ai,double((Nx-2)*(Ny-2)));
printf("Mean curvature = %f (analytical =0) \n", plane.Ji);
printf("Euler characteristic = %f (analytical = 2.0) \n",plane.Xi);
} }
MPI_Barrier(comm); MPI_Barrier(comm);
MPI_Finalize(); MPI_Finalize();