More debugging

Former-commit-id: b9e8ae9762c903ae2669c9c3562fb93b887180f4

include/igl/shapeup.cpp

@@ -68,6 +68,13 @@ namespace igl
             DCloseTriplets.push_back(Triplet<double>(i,b(i), 1.0));
         
         sudata.DClose.setFromTriplets(DCloseTriplets.begin(), DCloseTriplets.end());
+
+		//Building smoothness matrix
+		std::vector<Triplet<double> > DSmoothTriplets;
+		for (int i = 0; i < E.rows(); i++) {
+			DSmoothTriplets.push_back(Triplet<double>(i, E(i, 0), -1));
+			DSmoothTriplets.push_back(Triplet<double>(i, E(i, 1), 1));
+		}
         
         igl::cat(1, sudata.DShape, sudata.DClose, sudata.A);
         //is this allowed? repeating A.
@@ -95,7 +102,7 @@ namespace igl
         sudata.W.conservativeResize(SC.sum()+b.size()+E.rows(), SC.sum()+b.size()+E.rows());
         sudata.W.setFromTriplets(WTriplets.begin(), WTriplets.end());
         
-        sudata.At=sudata.A.transpose();
+        sudata.At=sudata.A.transpose();  //for efficieny, as we use the transpose a lot in the least squares
         sudata.Q=sudata.At*sudata.W*sudata.A;
         
         return min_quad_with_fixed_precompute(sudata.Q,VectorXi(),SparseMatrix<double>(),true,sudata.solver_data);
@@ -112,30 +119,31 @@ namespace igl
     {
         using namespace Eigen;
         using namespace std;
-        MatrixXd currP;
+        MatrixXd currP=P0;
         MatrixXd prevP=P0;
         MatrixXd projP;
-        MatrixXd b(sudata.A.rows(),3);
-        b.block(sudata.Q.rows(), 0, sudata.b.rows(),3)=bc;  //this stays constant throughout the iterations
+		MatrixXd rhs(sudata.A.rows(), 3); rhs.setZero();
+        rhs.block(sudata.DShape.rows(), 0, sudata.b.rows(),3)=bc;  //this stays constant throughout the iterations
         
         projP.conservativeResize(sudata.SC.rows(), 3*sudata.SC.maxCoeff());
         for (int i=0;i<sudata.maxIterations;i++){
             
             for (int j=0;j<sudata.SC.rows();j++)
                 sudata.local_projection(currP, sudata.SC,sudata.S,projP);
-            //constructing the projection part of the right hand side
+            //constructing the projection part of the (DShape rows of the) right hand side
             int currRow=0;
             for (int i=0;i<sudata.S.rows();i++){
                 for (int j=0;j<sudata.SC(i);j++){
-                    b.row(currRow++)=projP.block(i, 3*j, 1,3);
+					rhs.row(currRow++)=projP.block(i, 3*j, 1,3);
                 }
             }
    
             //the global solve is independent per dimension
-            Eigen::PlainObjectBase<DerivedP> rhs=-sudata.At*sudata.W*b;
+            Eigen::PlainObjectBase<DerivedP> lsrhs=-sudata.At*sudata.W*rhs;
             min_quad_with_fixed_solve(sudata.solver_data, rhs,Eigen::PlainObjectBase<DerivedP>(),Eigen::PlainObjectBase<DerivedP>(), currP);
 
             double currChange=(currP-prevP).lpNorm<Infinity>();
+			cout << "currChange: " << currChange << endl;
             prevP=currP;
             if (currChange<sudata.pTolerance)
                 break;

tutorial/801_ShapeUp/main.cpp

@@ -1,7 +1,7 @@
 #include <igl/avg_edge_length.h>
 #include <igl/barycenter.h>
 #include <igl/jet.h>
-#include <igl/shapeup.h>
+#include <igl/shapeup.h>
 #include <igl/shapeup_local_projections.h>
 #include <igl/quad_planarity.h>
 #include <igl/readDMAT.h>
@@ -15,16 +15,16 @@
 
 // Quad mesh loaded
 Eigen::MatrixXd VQC;
-Eigen::MatrixXi FQC;
+Eigen::MatrixXi FQC;
 Eigen::MatrixXi E;
 Eigen::MatrixXi FQCtri;
-Eigen::MatrixXd PQC0, PQC1, PQC2, PQC3;
-// Euclidean-regular quad mesh
-Eigen::MatrixXd VQCregular;
-Eigen::MatrixXi FQCtriregular;
-Eigen::MatrixXd PQC0regular, PQC1regular, PQC2regular, PQC3regular;
-
-igl::ShapeupData su_data;
+Eigen::MatrixXd PQC0, PQC1, PQC2, PQC3;
+// Euclidean-regular quad mesh
+Eigen::MatrixXd VQCregular;
+Eigen::MatrixXi FQCtriregular;
+Eigen::MatrixXd PQC0regular, PQC1regular, PQC2regular, PQC3regular;
+
+igl::ShapeupData su_data;
 
 
 
@@ -39,7 +39,7 @@ bool key_down(igl::viewer::Viewer& viewer, unsigned char key, int modifier)
 
   // Plot the original quad mesh
   if (key == '1')
-  {
+  {
       cout<<"before setting mesh 1"<<endl;
     // Draw the triangulated quad mesh
     viewer.data.set_mesh(VQC, FQCtri);
@@ -63,7 +63,7 @@ bool key_down(igl::viewer::Viewer& viewer, unsigned char key, int modifier)
   // Plot the planarized quad mesh
   if (key == '2')
   {
-    // Draw the triangulated quad mesh
+    // Draw the triangulated quad mesh
       cout<<"before setting mesh 2"<<endl;
     viewer.data.set_mesh(VQCregular, FQCtri);
 
@@ -92,7 +92,7 @@ int main(int argc, char *argv[])
   using namespace std;
 
   // Load a quad mesh
-  igl::readOFF(TUTORIAL_SHARED_PATH "/halftunnel.off", VQC, FQC);
+  igl::readOFF(TUTORIAL_SHARED_PATH "/halftunnel.off", VQC, FQC);
 
   // Convert it in a triangle mesh
   FQCtri.resize(2*FQC.rows(), 3);
@@ -104,24 +104,24 @@ int main(int argc, char *argv[])
   igl::slice( VQC, FQC.col(3).eval(), 1, PQC3);
 
   // Create a planar version with ShapeUp
-  //igl::planarize_quad_mesh(VQC, FQC, 100, 0.005, VQCregular);
-    
-  E.resize(FQC.size(),2);
-  E.col(0)<<FQC.col(0),FQC.col(1),FQC.col(2),FQC.col(3);
-  E.col(1)<<FQC.col(1),FQC.col(2),FQC.col(3),FQC.col(0);
-    
-  VectorXi b;
-  VectorXd w(FQC.rows());
-  MatrixXd bc;
-    
-  VectorXi array_of_fours=VectorXi::Constant(FQC.rows(),4);
-    cout<<"before pre-computation"<<endl;
-    std::function<bool(const Eigen::PlainObjectBase<MatrixXd>&, const Eigen::PlainObjectBase<VectorXi>&, const Eigen::PlainObjectBase<MatrixXi>&, Eigen::PlainObjectBase<MatrixXd>&)> localFunction(igl::shapeup_identity_projection);
-    
-    shapeup_precomputation(VQC, array_of_fours,FQC,E,b,w, localFunction,su_data);
-    cout<<"after pre-computation"<<endl;
-    shapeup_solve(bc,VQC,su_data,VQCregular);
-    cout<<"after computation"<<endl;
+  //igl::planarize_quad_mesh(VQC, FQC, 100, 0.005, VQCregular);
+    
+  E.resize(FQC.size(),2);
+  E.col(0)<<FQC.col(0),FQC.col(1),FQC.col(2),FQC.col(3);
+  E.col(1)<<FQC.col(1),FQC.col(2),FQC.col(3),FQC.col(0);
+    
+  VectorXi b;
+  VectorXd w(FQC.rows());
+  MatrixXd bc(0,3);
+    
+  VectorXi array_of_fours=VectorXi::Constant(FQC.rows(),4);
+    cout<<"before pre-computation"<<endl;
+    std::function<bool(const Eigen::PlainObjectBase<MatrixXd>&, const Eigen::PlainObjectBase<VectorXi>&, const Eigen::PlainObjectBase<MatrixXi>&, Eigen::PlainObjectBase<MatrixXd>&)> localFunction(igl::shapeup_identity_projection);
+    
+    shapeup_precomputation(VQC, array_of_fours,FQC,E,b,w, localFunction,su_data);
+    cout<<"after pre-computation"<<endl;
+    shapeup_solve(bc,VQC,su_data,VQCregular);
+    cout<<"after computation"<<endl;
     
 
   // Convert the planarized mesh to triangles