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

debugging decl

Browse Source
This commit is contained in:
James E McClure
2018-09-13 22:31:19 -04:00
parent 33eeda579a
commit 3c953c8ab2
3 changed files with 32 additions and 10 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
analysis/Minkowski.cpp
View File

@@ -207,14 +207,14 @@ void Minkowski::ComputeScalar(const DoubleArray Field, const double isovalue)
a2 = object.EdgeAngle(e2);
a3 = object.EdgeAngle(e3);
Ji += 0.08333333333333*(a1*s1+a2*s2+a3*s3);
if (0.08333333333333*(a1*s1+a2*s2+a3*s3) < 0.f){
//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(" P={%f,%f,%f}\n",P1.x,P1.y,P1.z);
printf(" Q={%f,%f,%f}\n",P2.x,P2.y,P2.z);
printf(" R={%f,%f,%f}\n",P3.x,P3.y,P3.z);
}
//}
// Euler characteristic (half edge rule: one face - 0.5*(three edges))
Xi -= 0.5;
}

11
analysis/decl.cpp
View File

@@ -365,6 +365,7 @@ double DECL::EdgeAngle(int edge)
V = TriNormal(halfedge.twin(edge));
double dotprod=U.x*V.x + U.y*V.y + U.z*V.z;
if (dotprod > 1.f) dotprod=1.f;
if (dotprod < -1.f) dotprod=-1.f;
angle = acos(dotprod);
if (halfedge.twin(edge) < 0 ){
// turn in outward normal at cube face is pi/2 from each side of the cube
@@ -388,11 +389,11 @@ double DECL::EdgeAngle(int edge)
angle *= 0.5; // half edge value
}
//1.570796326794897
/*if (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(" 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);
}
*/
// 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);
//}
return angle;
}

27
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", { 32, 32, 32 } );
db->putVector<int>( "n", { 16, 16, 16 } );
db->putScalar<int>( "nspheres", 1 );
db->putVector<double>( "L", { 1, 1, 1 } );
return db;
@@ -60,13 +60,34 @@ int main(int argc, char **argv)
printf("Mean curvature = %f (analytical =0) \n", plane.Ji);
printf("Euler characteristic = %f (analytical = 2.0) \n",plane.Xi);
Minkowski cylinder(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) = 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;
}
}
}
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);
/*
Minkowski sphere(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) = 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)+(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;
}
}
}
@@ -77,7 +98,7 @@ int main(int argc, char **argv)
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);
*/
}
MPI_Barrier(comm);
MPI_Finalize();