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

mean curvautre getting closer

Browse Source
This commit is contained in:
James E McClure 2018-09-14 16:04:35 -04:00
parent 80441c7b4f
commit bbc5aaeb0f
3 changed files with 30 additions and 38 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
analysis/Minkowski.cpp
View File

@ -208,12 +208,11 @@ void Minkowski::ComputeScalar(const DoubleArray Field, const double isovalue)
a3 = object.EdgeAngle(e3);
Ji += 0.08333333333333*(a1*s1+a2*s2+a3*s3);
//if (0.08333333333333*(a1*s1+a2*s2+a3*s3) < 0.f){
double intcurv=0.08333333333333*(a1*s1+a2*s2+a3*s3);
double surfarea=sqrt(s*(s-s1)*(s-s2)*(s-s3));
//printf("%f, %f: s1=%f,s2=%f,s3=%f, a1=%f,a2=%f,a3=%f\n",surfarea,intcurv,s1,s2,s3,a1,a2,a3);
printf(" (%i,%i,%i) PQ(%i,%i)={%f,%f,%f} {%f,%f,%f} a=%f l=%f \n",i,j,k,e1,object.halfedge.twin(e1),P1.x,P1.y,P1.z,P2.x,P2.y,P2.z,a1,s1);
printf(" (%i,%i,%i) QR(%i,%i)={%f,%f,%f} {%f,%f,%f} a=%f l=%f \n",i,j,k,e2,object.halfedge.twin(e2),P2.x,P2.y,P2.z,P3.x,P3.y,P3.z,a2,s2);
printf(" (%i,%i,%i) RP(%i,%i)={%f,%f,%f} {%f,%f,%f} a=%f l=%f \n",i,j,k,e3,object.halfedge.twin(e3),P3.x,P3.y,P3.z,P1.x,P1.y,P1.z,a3,s3);
//double intcurv=0.08333333333333*(a1*s1+a2*s2+a3*s3);
//double surfarea=sqrt(s*(s-s1)*(s-s2)*(s-s3));
//printf(" (%i,%i,%i) PQ(%i,%i)={%f,%f,%f} {%f,%f,%f} a=%f l=%f \n",i,j,k,e1,object.halfedge.twin(e1),P1.x,P1.y,P1.z,P2.x,P2.y,P2.z,a1,s1);
//printf(" (%i,%i,%i) QR(%i,%i)={%f,%f,%f} {%f,%f,%f} a=%f l=%f \n",i,j,k,e2,object.halfedge.twin(e2),P2.x,P2.y,P2.z,P3.x,P3.y,P3.z,a2,s2);
//printf(" (%i,%i,%i) RP(%i,%i)={%f,%f,%f} {%f,%f,%f} a=%f l=%f \n",i,j,k,e3,object.halfedge.twin(e3),P3.x,P3.y,P3.z,P1.x,P1.y,P1.z,a3,s3);
//}
// Euler characteristic (half edge rule: one face - 0.5*(three edges))
Xi -= 0.5;

7
analysis/decl.cpp
View File

@ -396,12 +396,7 @@ double DECL::EdgeAngle(int edge)
}
angle *= 0.5; // half edge value
}
//1.570796326794897
// if (angle < 0.f){
printf(" %f, %f (Edge=%i, twin=%i)\n 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(" %f, %f (Edge=%i, twin=%i)\n 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);
return angle;
}

46
tests/TestSphereCurvature.cpp
View File

@ -17,7 +17,7 @@ std::shared_ptr<Database> loadInputs( )
auto db = std::make_shared<Database>();
db->putScalar<int>( "BC", 0 );
db->putVector<int>( "nproc", { 1, 1, 1 } );
db->putVector<int>( "n", { 16, 16, 16 } );
db->putVector<int>( "n", { 32, 32, 32 } );
db->putScalar<int>( "nspheres", 1 );
db->putVector<double>( "L", { 1, 1, 1 } );
return db;
@ -27,8 +27,8 @@ int main(int argc, char **argv)
{
MPI_Init(&argc,&argv);
MPI_Comm comm = MPI_COMM_WORLD;
int rank = MPI_WORLD_RANK();
int nprocs = MPI_WORLD_SIZE();
//int rank = MPI_WORLD_RANK();
//int nprocs = MPI_WORLD_SIZE();
int toReturn = 0;
{
int i,j,k;
@ -42,9 +42,10 @@ int main(int argc, char **argv)
Nx+=2; Ny+=2; Nz+=2;
DoubleArray Phase(Nx,Ny,Nz);
Minkowski plane(Dm);
printf("Set distance map \n");
printf("Set distance map for plane \n");
for (k=0; k<Nz; k++){
for (j=0; j<Ny; j++){
for (i=0; i<Nx; i++){
@ -53,52 +54,49 @@ int main(int argc, char **argv)
}
}
printf("Construct local isosurface \n");
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);
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 = 0) \n",plane.Xi);
Minkowski cylinder(Dm);
printf("Set distance map \n");
printf("Set distance map for cylinder \n");
for (k=0; k<Nz; k++){
for (j=0; j<Ny; j++){
for (i=0; i<Nx; i++){
//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))-0.3*Nx;
Phase(i,j,k) = sqrt((1.0*i-0.5*Nx)*(1.0*i-0.5*Nx)+(1.0*k-0.5*Nz)*(1.0*k-0.5*Nz))-0.3*Nx;
Phase(i,j,k) = sqrt((1.0*i-0.5*Nx)*(1.0*i-0.5*Nx)+(1.0*k-0.5*Nz)*(1.0*k-0.5*Nz))-0.4*Nx;
}
}
}
printf("Construct local isosurface \n");
printf(" Construct local isosurface \n");
cylinder.ComputeScalar(Phase,0.f);
printf("Surface area = %f (analytical = %f) \n", cylinder.Ai,2*3.14159*0.3*double((Nx-2)*Nx));
printf("Mean curvature = %f (analytical = %f) \n", cylinder.Ji,2*3.14159*double((Nx-2)));
printf("Euler characteristic = %f (analytical = 2.0) \n",cylinder.Xi);
printf(" Surface area = %f (analytical = %f) \n", cylinder.Ai,2*3.14159*0.4*double((Nx-2)*Nx));
printf(" Mean curvature = %f (analytical = %f) \n", cylinder.Ji,2*3.14159*double((Nx-2)));
printf(" Euler characteristic = %f (analytical = 0) \n",cylinder.Xi);
/*
Minkowski sphere(Dm);
printf("Set distance map \n");
printf("Set distance map for a sphere \n");
for (k=0; k<Nz; k++){
for (j=0; j<Ny; j++){
for (i=0; i<Nx; i++){
// 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;
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))-0.3*Nx;
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.4*Nx;
}
}
}
printf("Construct local isosurface \n");
printf(" Construct local isosurface \n");
sphere.ComputeScalar(Phase,0.f);
printf("Surface area = %f (analytical = %f) \n", sphere.Ai,4*3.14159*0.3*0.3*double(Nx*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(" Surface area = %f (analytical = %f) \n", sphere.Ai,4*3.14159*0.16*double(Nx*Nx));
printf(" Mean curvature = %f (analytical = %f) \n", sphere.Ji,8*3.14159*0.4*double(Nx));
printf(" Euler characteristic = %f (analytical = 2.0) \n",sphere.Xi);
}
MPI_Barrier(comm);
MPI_Finalize();