From bbc5aaeb0f2a1e0a50261bf71b6a8f004790a518 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Fri, 14 Sep 2018 16:04:35 -0400 Subject: [PATCH] mean curvautre getting closer --- analysis/Minkowski.cpp | 11 ++++---- analysis/decl.cpp | 7 +---- tests/TestSphereCurvature.cpp | 50 +++++++++++++++++------------------ 3 files changed, 30 insertions(+), 38 deletions(-) diff --git a/analysis/Minkowski.cpp b/analysis/Minkowski.cpp index a2dd191a..61a9463c 100644 --- a/analysis/Minkowski.cpp +++ b/analysis/Minkowski.cpp @@ -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; diff --git a/analysis/decl.cpp b/analysis/decl.cpp index 95288e02..e3f19f32 100644 --- a/analysis/decl.cpp +++ b/analysis/decl.cpp @@ -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; } diff --git a/tests/TestSphereCurvature.cpp b/tests/TestSphereCurvature.cpp index 97e1c3c5..3f543c9f 100644 --- a/tests/TestSphereCurvature.cpp +++ b/tests/TestSphereCurvature.cpp @@ -17,7 +17,7 @@ std::shared_ptr loadInputs( ) auto db = std::make_shared(); db->putScalar( "BC", 0 ); db->putVector( "nproc", { 1, 1, 1 } ); - db->putVector( "n", { 16, 16, 16 } ); + db->putVector( "n", { 32, 32, 32 } ); db->putScalar( "nspheres", 1 ); db->putVector( "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