From 232db68c7bf45ecec428d83ce313ac45a2f321d3 Mon Sep 17 00:00:00 2001
From: Olga Diamanti <olga.diamanti@inf.ethz.ch>
Date: Thu, 16 Jul 2015 17:50:03 +0200
Subject: [PATCH] tutorial for integrable polyvectors

---
 include/igl/integrable_polyvector_fields.cpp  | 243 +++---
 include/igl/integrable_polyvector_fields.h    |   3 -
 ...tor_field_comb_from_matchings_and_cuts.cpp |   4 +-
 tutorial/705_Integrable/CMakeLists.txt        |  11 +
 tutorial/705_Integrable/main.cpp              | 709 ++++++++++++++++++
 tutorial/CMakeLists.txt                       |   1 +
 6 files changed, 843 insertions(+), 128 deletions(-)
 create mode 100644 tutorial/705_Integrable/CMakeLists.txt
 create mode 100644 tutorial/705_Integrable/main.cpp

diff --git a/include/igl/integrable_polyvector_fields.cpp b/include/igl/integrable_polyvector_fields.cpp
index cec1bba84..be356829a 100644
--- a/include/igl/integrable_polyvector_fields.cpp
+++ b/include/igl/integrable_polyvector_fields.cpp
@@ -29,8 +29,7 @@ wCloseUnconstrained(1e-3),
 wCloseConstrained(100),
 redFactor_wsmooth(.8),
 gamma(0.1),
-tikh_gamma(1e-8),
-do_partial(false)
+tikh_gamma(1e-8)
 {}
 
 
@@ -41,7 +40,7 @@ namespace igl {
   class IntegrableFieldSolver
   {
   private:
-    
+
     IntegrableFieldSolverData<DerivedV, DerivedF, DerivedFF, DerivedC> &data;
     //Symbolic calculations
     IGL_INLINE void rj_barrier_face(const Eigen::RowVectorXd &vec2D_a,
@@ -71,24 +70,24 @@ namespace igl {
                                        Eigen::VectorXd &residuals,
                                        bool do_jac = false,
                                        Eigen::MatrixXd &Jac = *(Eigen::MatrixXd*)NULL);
-    
+
   public:
     IGL_INLINE IntegrableFieldSolver(IntegrableFieldSolverData<DerivedV, DerivedF, DerivedFF, DerivedC> &cffsoldata);
-    
+
     IGL_INLINE bool solve(integrable_polyvector_fields_parameters &params,
                           Eigen::PlainObjectBase<DerivedFF>& current_field,
                           bool current_field_is_not_ccw);
-    
+
     IGL_INLINE void solveGaussNewton(integrable_polyvector_fields_parameters &params,
                                      const Eigen::VectorXd &x_initial,
                                      Eigen::VectorXd &x);
-    
+
     //Compute residuals and Jacobian for Gauss Newton
     IGL_INLINE double RJ(const Eigen::VectorXd &x,
                          const Eigen::VectorXd &x0,
                          const integrable_polyvector_fields_parameters &params,
                          bool doJacs = false);
-    
+
     IGL_INLINE void RJ_Smoothness(const Eigen::MatrixXd &sol2D,
                                   const double &wSmoothSqrt,
                                   const int startRowInJacobian,
@@ -104,7 +103,6 @@ namespace igl {
                                  const Eigen::MatrixXd &sol02D,
                                  const double &wCloseUnconstrainedSqrt,
                                  const double &wCloseConstrainedSqrt,
-                                 const bool do_partial,
                                  const int startRowInJacobian,
                                  bool doJacs = false,
                                  const int startIndexInVectors = 0);
@@ -119,8 +117,8 @@ namespace igl {
                                 const int startRowInJacobian,
                                 bool doJacs = false,
                                 const int startIndexInVectors = 0);
-    
-    
+
+
   };
 };
 
@@ -177,7 +175,7 @@ initializeConstraints(const Eigen::VectorXi& b,
   //sort in descending order (so removing rows will work)
   igl::sort(constrained_unsorted, 1, false, constrained, iSorted);
   constrained_vec3 = bc.template cast<double>();
-  
+
 }
 
 template<typename DerivedV, typename DerivedF, typename DerivedFF, typename DerivedC>
@@ -214,7 +212,7 @@ makeFieldCCW(Eigen::MatrixXd &sol3D)
       int n1 = 3*2*2-n2;
       t.segment(0,n1) = t.segment(s1,n1);
       t.segment(n1, n2) = temp;
-      
+
     }
     sol3D.row(fi) = t.segment(0, 2*3);
   }
@@ -231,7 +229,7 @@ initializeOriginalVariable(const Eigen::PlainObjectBase<DerivedFF>& original_fie
   xOriginal.setZero(numVariables);
   for (int i =0; i<numF; i++)
     xOriginal.segment(i*2*2, 2*2) = sol2D.row(i);
-  
+
 }
 
 template<typename DerivedV, typename DerivedF, typename DerivedFF, typename DerivedC>
@@ -243,7 +241,7 @@ computeInteriorEdges()
   numInteriorEdges = 0;
   isBorderEdge.setZero(numE,1);
   indFullToInterior = -1.*Eigen::VectorXi::Ones(numE,1);
-  
+
   for(unsigned i=0; i<numE; ++i)
   {
     if ((E2F(i,0) == -1) || ((E2F(i,1) == -1)))
@@ -254,7 +252,7 @@ computeInteriorEdges()
       numInteriorEdges++;
     }
   }
-  
+
   E2F_int.resize(numInteriorEdges, 2);
   indInteriorToFull.setZero(numInteriorEdges,1);
   int ii = 0;
@@ -266,7 +264,7 @@ computeInteriorEdges()
     indInteriorToFull[ii] = k;
     ii++;
   }
-  
+
 }
 
 template<typename DerivedV, typename DerivedF, typename DerivedFF, typename DerivedC>
@@ -294,7 +292,7 @@ add_jac_indices_face(const int numInnerRows,
 {
   for (int fi=0; fi<numF; ++fi)
   {
-    
+
     int startRow = startRowInJacobian+numInnerRows*fi;
     int startIndex = startIndexInVectors+numInnerRows*numInnerCols*fi;
     face_Jacobian_indices(startRow, startIndex, fi, 2, numInnerRows, numInnerCols, Rows, Columns);
@@ -340,13 +338,13 @@ add_jac_indices_edge(const int numInnerRows,
     // the two faces of the flap
     int a = E2F_int(ii,0);
     int b = E2F_int(ii,1);
-    
-    
+
+
     int startRow = startRowInJacobian+numInnerRows*ii;
     int startIndex = startIndexInVectors+numInnerRows*numInnerCols*ii;
-    
+
     edge_Jacobian_indices(startRow, startIndex, a, b, 2, numInnerRows, numInnerCols, Rows, Columns);
-    
+
   }
 }
 
@@ -386,30 +384,30 @@ computeJacobianPattern()
   num_residuals_polycurl = 2*numInteriorEdges;
   num_residuals_quotcurl = numInteriorEdges;
   num_residuals_barrier = numF;
-  
+
   num_residuals = num_residuals_smooth + num_residuals_close + num_residuals_polycurl + num_residuals_quotcurl + num_residuals_barrier;
-  
+
   residuals.setZero(num_residuals,1);
-  
+
   //per edge
   numJacElements_smooth = numInteriorEdges*numInnerJacCols_edge*numInnerJacRows_smooth;
   numJacElements_polycurl = numInteriorEdges*numInnerJacCols_edge*numInnerJacRows_polycurl;
   numJacElements_quotcurl = numInteriorEdges*numInnerJacCols_edge*numInnerJacRows_quotcurl;
-  
+
   //per face
   numJacElements_barrier = numF*numInnerJacCols_face*numInnerJacRows_barrier;
   numJacElements_close = numF*numInnerJacCols_face*numInnerJacRows_close;
-  
+
   numJacElements = (numJacElements_smooth +numJacElements_polycurl + numJacElements_quotcurl) + (numJacElements_barrier +numJacElements_close);
   //allocate
   II_Jac.setZero(numJacElements);
   JJ_Jac.setZero(numJacElements);
   SS_Jac.setOnes(numJacElements);
-  
+
   //set stuff (attention: order !)
   int startRowInJacobian = 0;
   int startIndexInVectors = 0;
-  
+
   //smoothness
   add_jac_indices_edge(numInnerJacRows_smooth,
                        numInnerJacCols_edge,
@@ -419,7 +417,7 @@ computeJacobianPattern()
                        JJ_Jac);
   startRowInJacobian += num_residuals_smooth;
   startIndexInVectors += numJacElements_smooth;
-  
+
   //closeness
   add_jac_indices_face(numInnerJacRows_close,
                        numInnerJacCols_face,
@@ -429,7 +427,7 @@ computeJacobianPattern()
                        JJ_Jac);
   startRowInJacobian += num_residuals_close;
   startIndexInVectors += numJacElements_close;
-  
+
   //barrier
   add_jac_indices_face(numInnerJacRows_barrier,
                        numInnerJacCols_face,
@@ -439,7 +437,7 @@ computeJacobianPattern()
                        JJ_Jac);
   startRowInJacobian += num_residuals_barrier;
   startIndexInVectors += numJacElements_barrier;
-  
+
   //polycurl
   add_jac_indices_edge(numInnerJacRows_polycurl,
                        numInnerJacCols_edge,
@@ -449,7 +447,7 @@ computeJacobianPattern()
                        JJ_Jac);
   startRowInJacobian += num_residuals_polycurl;
   startIndexInVectors += numJacElements_polycurl;
-  
+
   //quotcurl
   add_jac_indices_edge(numInnerJacRows_quotcurl,
                        numInnerJacCols_edge,
@@ -505,7 +503,7 @@ computeNewHessValues()
                                               Hess_triplets[i].col(),
                                               SS_Jac(indInSS_Hess_1_vec[i])*SS_Jac(indInSS_Hess_2_vec[i])
                                               );
-  
+
   Hess.setFromTriplets(Hess_triplets.begin(), Hess_triplets.end());
 }
 
@@ -526,13 +524,13 @@ solve(igl::integrable_polyvector_fields_parameters &params,
   Eigen::MatrixXd sol3D = current_field.template cast<double>();
   if (current_field_is_not_ccw)
     data.makeFieldCCW(sol3D);
-  
+
   igl::global2local(data.B1, data.B2, sol3D, sol2D);
   Eigen::VectorXd x;
   x.setZero(data.numVariables);
   for (int i =0; i<data.numF; i++)
     x.segment(i*2*2, 2*2) = sol2D.row(i);
-  
+
   //get x
   solveGaussNewton(params, data.xOriginal, x);
   //get output from x
@@ -550,23 +548,23 @@ solveGaussNewton(integrable_polyvector_fields_parameters &params,
                  Eigen::VectorXd &x)
 {
   bool converged = false;
-  
+
   double F;
   Eigen::VectorXd xprev = x;
   Eigen::VectorXd xc = igl::slice(x_initial, data.constrained, 1);
   //  double ESmooth, EClose, ECurl, EQuotCurl, EBarrier;
   for (int innerIter = 0; innerIter<params.numIter; ++innerIter)
   {
-    
+
     //set constrained entries to those of the initial
     igl::slice_into(xc, data.constrained, 1, xprev);
-    
-    
+
+
     //get function, gradients and Hessians
     F = RJ(x, xprev, params, true);
-    
+
     printf("IntegrableFieldSolver -- Iteration %d\n", innerIter);
-    
+
     if((data.residuals.array() == std::numeric_limits<double>::infinity()).any())
     {
       std::cerr<<"IntegrableFieldSolver -- residuals: got infinity somewhere"<<std::endl;
@@ -577,20 +575,20 @@ solveGaussNewton(integrable_polyvector_fields_parameters &params,
       std::cerr<<"IntegrableFieldSolver -- residuals: got infinity somewhere"<<std::endl;
       exit(1);
     };
-    
+
     converged = false;
-    
+
     Eigen::VectorXd rhs = data.Jac.transpose()*data.residuals;
-    
+
     bool success;
     data.solver.factorize(data.Hess);
     success = data.solver.info() == Eigen::Success;
-    
+
     if(!success)
       std::cerr<<"IntegrableFieldSolver -- Could not do LU"<<std::endl;
-    
+
     Eigen::VectorXd direction;
-    
+
     double error;
     direction = data.solver.solve(rhs);
     error = (data.Hess*direction - rhs).cwiseAbs().maxCoeff();
@@ -598,7 +596,7 @@ solveGaussNewton(integrable_polyvector_fields_parameters &params,
     {
       std::cerr<<"IntegrableFieldSolver -- Could not solve"<<std::endl;
     }
-    
+
     // adaptive backtracking
     bool repeat = true;
     int run = 0;
@@ -626,8 +624,8 @@ solveGaussNewton(integrable_polyvector_fields_parameters &params,
       }
       run++;
     }
-    
-    
+
+
     if (!converged)
     {
       xprev = x;
@@ -639,8 +637,8 @@ solveGaussNewton(integrable_polyvector_fields_parameters &params,
       break;
     }
   }
-  
-  
+
+
 }
 
 template <typename DerivedV, typename DerivedF, typename DerivedFF, typename DerivedC>
@@ -655,43 +653,43 @@ RJ(const Eigen::VectorXd &x,
     sol2D.row(i) = x.segment(i*2*2, 2*2);
   for (int i =0; i<data.numF; i++)
     sol02D.row(i) = x0.segment(i*2*2, 2*2);
-  
+
   data.SS_Jac.setZero(data.numJacElements);
-  
+
   //set stuff (attention: order !)
   int startRowInJacobian = 0;
   int startIndexInVectors = 0;
-  
+
   //smoothness
   RJ_Smoothness(sol2D, sqrt(params.wSmooth), startRowInJacobian, doJacs, startIndexInVectors);
   startRowInJacobian += data.num_residuals_smooth;
   startIndexInVectors += data.numJacElements_smooth;
-  
+
   //closeness
-  RJ_Closeness(sol2D, sol02D, sqrt(params.wCloseUnconstrained), sqrt(params.wCloseConstrained), params.do_partial, startRowInJacobian, doJacs, startIndexInVectors);
+  RJ_Closeness(sol2D, sol02D, sqrt(params.wCloseUnconstrained), sqrt(params.wCloseConstrained), startRowInJacobian, doJacs, startIndexInVectors);
   startRowInJacobian += data.num_residuals_close;
   startIndexInVectors += data.numJacElements_close;
-  
+
   //barrier
   RJ_Barrier(sol2D, params.sBarrier, sqrt(params.wBarrier), startRowInJacobian, doJacs, startIndexInVectors);
   startRowInJacobian += data.num_residuals_barrier;
   startIndexInVectors += data.numJacElements_barrier;
-  
+
   //polycurl
   RJ_Curl(sol2D, params.wCurl, powl(params.wCurl, 2), startRowInJacobian, doJacs, startIndexInVectors);
   startRowInJacobian += data.num_residuals_polycurl;
   startIndexInVectors += data.numJacElements_polycurl;
-  
+
   //quotcurl
   RJ_QuotCurl(sol2D, sqrt(params.wQuotCurl), startRowInJacobian, doJacs, startIndexInVectors);
-  
+
   if(doJacs)
   {
     for (int i =0; i<data.numJacElements; ++i)
       data.Jac.coeffRef(data.II_Jac(i), data.JJ_Jac(i)) = data.SS_Jac(i);
     data.computeNewHessValues();
   }
-  
+
   return data.residuals.transpose()*data.residuals;
 }
 
@@ -714,7 +712,7 @@ rj_smoothness_edge(const Eigen::RowVectorXd &vec2D_a,
   const Eigen::RowVector2d &vb = vec2D_b.segment(2, 2);
   const double &xua=ua[0], &yua=ua[1], &xva=va[0], &yva=va[1];
   const double &xub=ub[0], &yub=ub[1], &xvb=vb[0], &yvb=vb[1];
-  
+
   double xua_2 = xua*xua;
   double xva_2 = xva*xva;
   double yua_2 = yua*yua;
@@ -723,12 +721,12 @@ rj_smoothness_edge(const Eigen::RowVectorXd &vec2D_a,
   double xvb_2 = xvb*xvb;
   double yub_2 = yub*yub;
   double yvb_2 = yvb*yvb;
-  
+
   double sA = sin(nA*k);
   double cA = cos(nA*k);
   double sB = sin(nB*k);
   double cB = cos(nB*k);
-  
+
   double t1 = xua*yua;
   double t2 = xva*yva;
   double t3 = xub*xvb;
@@ -737,33 +735,33 @@ rj_smoothness_edge(const Eigen::RowVectorXd &vec2D_a,
   double t6 = xub*yub;
   double t7 = yua*yva;
   double t8 = xvb*yvb;
-  
+
   double t9  = xva_2 - yva_2;
   double t10 = xua_2 - yua_2;
   double t11 = xvb_2 - yvb_2;
   double t12 = xub_2 - yub_2;
-  
+
   double t13 = 2*t1 + 2*t2;
-  
+
   double t17 = (2*t1*t9 + 2*t2*t10);
   double t19 = (t10*t9 - 4*t5*t7);
-  
-  
+
+
   residuals.resize(4, 1);
   residuals <<
   cA*(t10 + t9) - sA*(t13) - t12 - t11,
   sA*(t10 + t9) - 2*t8 - 2*t6 + cA*(t13),
   cB*t19 - (t12)*(t11) - sB*t17 + 4*t3*t4,
   cB*t17 + sB*t19 - 2*t6*(t11) - 2*t8*(t12);
-  
-  
+
+
   if (do_jac)
   {
     double t20 = 2*yua*t9 + 4*xua*t2;
     double t21 = 2*xua*t9 - 4*xva*t7;
     double t22 = 2*yva*t10 + 4*t5*yua;
     double t23 = 2*xva*t10 - 4*t1*yva;
-    
+
     Jac.resize(4,8);
     Jac <<                                                                     2*xua*cA - 2*yua*sA,                                                                     - 2*yua*cA - 2*xua*sA,                                                                     2*xva*cA - 2*yva*sA,                                                                     - 2*yva*cA - 2*xva*sA,                                  -2*xub,                                 2*yub,                                  -2*xvb,                                 2*yvb,
     2*yua*cA + 2*xua*sA,                                                                       2*xua*cA - 2*yua*sA,                                                                     2*yva*cA + 2*xva*sA,                                                                       2*xva*cA - 2*yva*sA,                                  -2*yub,                                -2*xub,                                  -2*yvb,                                -2*xvb,
@@ -788,9 +786,9 @@ RJ_Smoothness(const Eigen::MatrixXd &sol2D,
     // the two faces of the flap
     int a = data.E2F_int(ii,0);
     int b = data.E2F_int(ii,1);
-    
+
     int k = data.indInteriorToFull[ii];
-    
+
     Eigen::MatrixXd tJac;
     Eigen::VectorXd tRes;
     rj_smoothness_edge(sol2D.row(a),
@@ -801,10 +799,10 @@ RJ_Smoothness(const Eigen::MatrixXd &sol2D,
                        tRes,
                        doJacs,
                        tJac);
-    
+
     int startRow = startRowInJacobian+data.numInnerJacRows_smooth*ii;
     data.residuals.segment(startRow,data.numInnerJacRows_smooth) = wSmoothSqrt*tRes;
-    
+
     if(doJacs)
     {
       int startIndex = startIndexInVectors+data.numInnerJacRows_smooth*data.numInnerJacCols_edge*ii;
@@ -822,26 +820,26 @@ rj_barrier_face(const Eigen::RowVectorXd &vec2D_a,
                 bool do_jac,
                 Eigen::MatrixXd &Jac)
 {
-  
+
   const Eigen::RowVector2d &ua = vec2D_a.segment(0, 2);
   const Eigen::RowVector2d &va = vec2D_a.segment(2, 2);
   const double &xua=ua[0], &yua=ua[1], &xva=va[0], &yva=va[1];
-  
-  
+
+
   double xva_2 = xva*xva;
   double yua_2 = yua*yua;
   double xua_2 = xua*xua;
   double yva_2 = yva*yva;
-  
+
   double s_2 = s*s;
   double s_3 = s*s_2;
-  
+
   double t00 = xua*yva;
   double t01 = xva*yua;
   double t05 = t00 - t01;
   double t05_2 = t05*t05;
   double t05_3 = t05*t05_2;
-  
+
   if (do_jac)
     Jac.setZero(1,4);
   residuals.resize(1, 1);
@@ -875,7 +873,7 @@ RJ_Barrier(const Eigen::MatrixXd &sol2D,
 {
   if (wBarrierSqrt ==0)
     return;
-  
+
   for (int fi=0; fi<data.numF; ++fi)
   {
     Eigen::MatrixXd tJac;
@@ -885,10 +883,10 @@ RJ_Barrier(const Eigen::MatrixXd &sol2D,
                     tRes,
                     doJacs,
                     tJac);
-    
+
     int startRow = startRowInJacobian+ data.numInnerJacRows_barrier * fi;
     data.residuals.segment(startRow,data.numInnerJacRows_barrier) = wBarrierSqrt*tRes;
-    
+
     if(doJacs)
     {
       int startIndex = startIndexInVectors+data.numInnerJacRows_barrier*data.numInnerJacCols_face*fi;
@@ -903,7 +901,6 @@ RJ_Closeness(const Eigen::MatrixXd &sol2D,
              const Eigen::MatrixXd &sol02D,
              const double &wCloseUnconstrainedSqrt,
              const double &wCloseConstrainedSqrt,
-             const bool do_partial,
              const int startRowInJacobian,
              bool doJacs,
              const int startIndexInVectors)
@@ -917,27 +914,27 @@ RJ_Closeness(const Eigen::MatrixXd &sol2D,
       weights.setConstant(wCloseUnconstrainedSqrt);
     else
     {
-      if (do_partial&&data.is_constrained_face[fi]==1)
+      if (data.is_constrained_face[fi]==1)
         //only constrain the first vector
         weights<<wCloseConstrainedSqrt,wCloseConstrainedSqrt,wCloseUnconstrainedSqrt,wCloseUnconstrainedSqrt;
       else
         //either not partial, or partial with 2 constraints
         weights.setConstant(wCloseConstrainedSqrt);
     }
-    
+
     Eigen::MatrixXd tJac;
     Eigen::VectorXd tRes;
     tJac = Eigen::MatrixXd::Identity(4,4);
     tRes.resize(4, 1); tRes<<(sol2D.row(fi)-sol02D.row(fi)).transpose();
     int startRow = startRowInJacobian+data.numInnerJacRows_close*fi;
     data.residuals.segment(startRow,data.numInnerJacRows_close) = weights.array()*tRes.array();
-    
+
     if(doJacs)
     {
       int startIndex = startIndexInVectors+data.numInnerJacRows_close*data.numInnerJacCols_face*fi;
       data.add_Jacobian_to_svector(startIndex, weights.asDiagonal()*tJac,data.SS_Jac);
     }
-    
+
   }
 }
 
@@ -960,26 +957,26 @@ rj_polycurl_edge(const Eigen::RowVectorXd &vec2D_a,
   const double &xub=ub[0], &yub=ub[1], &xvb=vb[0], &yvb=vb[1];
   const double &xea=ea[0], &yea=ea[1];
   const double &xeb=eb[0], &yeb=eb[1];
-  
+
   const double dua = (xea*xua + yea*yua);
   const double dub = (xeb*xub + yeb*yub);
   const double dva = (xea*xva + yea*yva);
   const double dvb = (xeb*xvb + yeb*yvb);
-  
+
   const double dua_2 = dua*dua;
   const double dva_2 = dva*dva;
   const double dub_2 = dub*dub;
   const double dvb_2 = dvb*dvb;
-  
+
   residuals.resize(2, 1);
   residuals << dua_2 - dub_2 + dva_2 - dvb_2,
   dua_2*dva_2 - dub_2*dvb_2 ;
-  
-  
-  
+
+
+
   if (do_jac)
   {
-    
+
     Jac.resize(2,8);
     Jac << 2*xea*dua,                       2*yea*dua,                       2*xea*dva,                       2*yea*dva,                       -2*xeb*dub,                       -2*yeb*dub,                       -2*xeb*dvb,                       -2*yeb*dvb,
     2*xea*dua*dva_2, 2*yea*dua*dva_2, 2*xea*dua_2*dva, 2*yea*dua_2*dva, -2*xeb*dub*dvb_2, -2*yeb*dub*dvb_2, -2*xeb*dub_2*dvb, -2*yeb*dub_2*dvb;
@@ -1002,12 +999,12 @@ RJ_Curl(const Eigen::MatrixXd &sol2D,
     // the two faces of the flap
     int a = data.E2F_int(ii,0);
     int b = data.E2F_int(ii,1);
-    
+
     int k = data.indInteriorToFull[ii];
-    
+
     // the common edge, a 3 vector
     const Eigen::RowVector3d &ce = data.EVecNorm.row(k);
-    
+
     // the common edge expressed in local coordinates in the two faces
     // x/y denotes real/imaginary
     double xea = data.B1.row(a).dot(ce);
@@ -1016,8 +1013,8 @@ RJ_Curl(const Eigen::MatrixXd &sol2D,
     double xeb = data.B1.row(b).dot(ce);
     double yeb = data.B2.row(b).dot(ce);
     Eigen::RowVector2d eb; eb<<xeb, yeb;
-    
-    
+
+
     Eigen::MatrixXd tJac;
     Eigen::VectorXd tRes;
     rj_polycurl_edge(sol2D.row(a),
@@ -1027,14 +1024,14 @@ RJ_Curl(const Eigen::MatrixXd &sol2D,
                      tRes,
                      doJacs,
                      tJac);
-    
+
     tRes[0] = tRes[0]*wCASqrt;
     tRes[1] = tRes[1]*wCBSqrt;
-    
-    
+
+
     int startRow = startRowInJacobian+data.numInnerJacRows_polycurl*ii;
     data.residuals.segment(startRow,data.numInnerJacRows_polycurl) = tRes;
-    
+
     if(doJacs)
     {
       tJac.row(0) = tJac.row(0)*wCASqrt;
@@ -1042,7 +1039,7 @@ RJ_Curl(const Eigen::MatrixXd &sol2D,
       int startIndex = startIndexInVectors+data.numInnerJacRows_polycurl*data.numInnerJacCols_edge*ii;
       data.add_Jacobian_to_svector(startIndex, tJac,data.SS_Jac);
     }
-    
+
   }
 }
 
@@ -1065,23 +1062,23 @@ rj_quotcurl_edge_polyversion(const Eigen::RowVectorXd &vec2D_a,
   const double &xub=ub[0], &yub=ub[1], &xvb=vb[0], &yvb=vb[1];
   const double &xea=ea[0], &yea=ea[1];
   const double &xeb=eb[0], &yeb=eb[1];
-  
+
   double dua = (xea*xua + yea*yua);
   double dva = (xea*xva + yea*yva);
   double dub = (xeb*xub + yeb*yub);
   double dvb = (xeb*xvb + yeb*yvb);
-  
+
   double dua_2 = dua * dua;
   double dva_2 = dva * dva;
   double dub_2 = dub * dub;
   double dvb_2 = dvb * dvb;
   double t00 = (dub_2 - dvb_2);
   double t01 = (dua_2 - dva_2);
-  
-  
+
+
   residuals.resize(1, 1);
   residuals << dua*dva*t00 - dub*dvb*t01;
-  
+
   if (do_jac)
   {
     Jac.resize(1,8);
@@ -1102,12 +1099,12 @@ RJ_QuotCurl(const Eigen::MatrixXd &sol2D,
     // the two faces of the flap
     int a = data.E2F_int(ii,0);
     int b = data.E2F_int(ii,1);
-    
+
     int k = data.indInteriorToFull[ii];
-    
+
     // the common edge, a 3 vector
     const Eigen::RowVector3d &ce = data.EVecNorm.row(k);
-    
+
     // the common edge expressed in local coordinates in the two faces
     // x/y denotes real/imaginary
     double xea = data.B1.row(a).dot(ce);
@@ -1116,8 +1113,8 @@ RJ_QuotCurl(const Eigen::MatrixXd &sol2D,
     double xeb = data.B1.row(b).dot(ce);
     double yeb = data.B2.row(b).dot(ce);
     Eigen::RowVector2d eb; eb<<xeb, yeb;
-    
-    
+
+
     Eigen::MatrixXd tJac;
     Eigen::VectorXd tRes;
     rj_quotcurl_edge_polyversion(sol2D.row(a),
@@ -1127,10 +1124,10 @@ RJ_QuotCurl(const Eigen::MatrixXd &sol2D,
                                  tRes,
                                  doJacs,
                                  tJac);
-    
+
     int startRow = startRowInJacobian+ data.numInnerJacRows_quotcurl*ii;
     data.residuals.segment(startRow,data.numInnerJacRows_quotcurl) = wQuotCurlSqrt*tRes;
-    
+
     if(doJacs)
     {
       int startIndex = startIndexInVectors+data.numInnerJacRows_quotcurl*data.numInnerJacCols_edge*ii;
@@ -1151,11 +1148,11 @@ IGL_INLINE void igl::integrable_polyvector_fields_precompute(
                                                              igl::IntegrableFieldSolverData<DerivedV, DerivedF, DerivedFF, DerivedC> &data)
 {
   data.precomputeMesh(V,F);
-  
+
   data.computeJacobianPattern();
   data.computeHessianPattern();
   data.solver.analyzePattern(data.Hess);
-  
+
   data.initializeConstraints(b,bc,constraint_level);
   data.initializeOriginalVariable(original_field);
 };
diff --git a/include/igl/integrable_polyvector_fields.h b/include/igl/integrable_polyvector_fields.h
index 363744f3a..8869e1e35 100644
--- a/include/igl/integrable_polyvector_fields.h
+++ b/include/igl/integrable_polyvector_fields.h
@@ -101,9 +101,6 @@ struct igl::integrable_polyvector_fields_parameters
   double gamma;
   //tikhonov regularization term (typically not needed, default value should suffice)
   double tikh_gamma;
-  //boolean, determines whether the frames at the constrained faces should be treated
-  //as partial constraints (i.e. only one of the vectors is preserved in the result)
-  bool do_partial;
 
   IGL_INLINE integrable_polyvector_fields_parameters();
 
diff --git a/include/igl/polyvector_field_comb_from_matchings_and_cuts.cpp b/include/igl/polyvector_field_comb_from_matchings_and_cuts.cpp
index 1e8c5d09e..719bbcc2f 100644
--- a/include/igl/polyvector_field_comb_from_matchings_and_cuts.cpp
+++ b/include/igl/polyvector_field_comb_from_matchings_and_cuts.cpp
@@ -61,8 +61,8 @@ IGL_INLINE void igl::polyvector_field_comb_from_matchings_and_cuts(
       }
       else
       {
-        assert((E2Fcut(current_edge,1) == f0) &&
-               (E2Fcut(current_edge,0) == f1));
+        assert((E2F(current_edge,1) == f0) &&
+               (E2F(current_edge,0) == f1));
         //look at match_ba
         for(int i=0; i<half_degree; ++i)
         {
diff --git a/tutorial/705_Integrable/CMakeLists.txt b/tutorial/705_Integrable/CMakeLists.txt
new file mode 100644
index 000000000..63f380431
--- /dev/null
+++ b/tutorial/705_Integrable/CMakeLists.txt
@@ -0,0 +1,11 @@
+cmake_minimum_required(VERSION 2.6)
+project(705_Integrable)
+
+include("../CMakeLists.shared")
+
+set(SOURCES
+${PROJECT_SOURCE_DIR}/main.cpp
+)
+
+add_executable(${PROJECT_NAME}_bin ${SOURCES} ${SHARED_SOURCES})
+target_link_libraries(${PROJECT_NAME}_bin ${SHARED_LIBRARIES})
diff --git a/tutorial/705_Integrable/main.cpp b/tutorial/705_Integrable/main.cpp
new file mode 100644
index 000000000..0222dba82
--- /dev/null
+++ b/tutorial/705_Integrable/main.cpp
@@ -0,0 +1,709 @@
+#include <igl/readOFF.h>
+#include <igl/readOBJ.h>
+#include <igl/n_polyvector.h>
+#include <igl/integrable_polyvector_fields.h>
+#include <igl/viewer/Viewer.h>
+#include <igl/local_basis.h>
+#include <igl/avg_edge_length.h>
+#include <igl/is_border_vertex.h>
+#include <igl/adjacency_list.h>
+#include <igl/vertex_triangle_adjacency.h>
+#include <igl/triangle_triangle_adjacency.h>
+#include <igl/edge_topology.h>
+#include <igl/jet.h>
+#include <igl/barycenter.h>
+#include <igl/polyvector_field_matchings.h>
+#include <igl/polyvector_field_singularities_from_matchings.h>
+#include <igl/polyvector_field_cut_mesh_with_singularities.h>
+#include <igl/polyvector_field_comb_from_matchings_and_cuts.h>
+#include <igl/polyvector_field_poisson_reconstruction.h>
+#include <igl/cut_mesh.h>
+#include <igl/slice.h>
+#include <igl/false_barycentric_subdivision.h>
+
+#include <iostream>
+#include <fstream>
+#include <igl/matlab_format.h>
+using namespace std;
+
+// Input mesh
+Eigen::MatrixXd V;
+Eigen::MatrixXi F;
+std::vector<bool> V_border;
+std::vector<std::vector<int> > VF, VFi;
+std::vector<std::vector<int> > VV;
+Eigen::MatrixXi TT, TTi;
+Eigen::MatrixXi E, E2F, F2E;
+
+// Per face bases (only needed to generate constraints)
+Eigen::MatrixXd B1,B2,B3;
+
+// "Subdivided" mesh obtained by splitting each triangle into 3 (only needed for display)
+Eigen::MatrixXd Vbs;
+Eigen::MatrixXi Fbs;
+
+// Scale for visualizing the fields
+double global_scale;
+
+// Scale for visualizing textures
+double uv_scale;
+
+// Data for original PolyVector field
+Eigen::MatrixXd two_pv_ori; // field
+Eigen::VectorXi singularities_ori; // singularities
+Eigen::VectorXd curl_ori; // curl per edge
+Eigen::MatrixXi cuts_ori; // cut edges
+Eigen::MatrixXd two_pv_poisson_ori; // field after poisson integration
+Eigen::VectorXf poisson_error_ori; // poisson integration error
+Eigen::MatrixXd scalars_ori;
+Eigen::MatrixXd Vcut_ori;
+Eigen::MatrixXi Fcut_ori;
+
+// Data for curl-free PolyVector field
+Eigen::MatrixXd two_pv; // field
+Eigen::VectorXi singularities; // singularities
+Eigen::VectorXd curl; // curl per edge
+Eigen::MatrixXi cuts; // cut edges
+Eigen::MatrixXd two_pv_poisson; // field after poisson integration
+Eigen::VectorXf poisson_error; // poisson integration error
+Eigen::MatrixXd scalars;
+Eigen::MatrixXd Vcut;
+Eigen::MatrixXi Fcut;
+
+// Vector of constrained faces
+Eigen::VectorXi b;
+
+// Matrix of constraints
+Eigen::MatrixXd bc;
+
+// "constraint level" flag (level=2 indicates that both directions are constrained,
+// level = 1 indicates a partially constrained face, i.e. only the first vector will
+// be constrained)
+Eigen::VectorXi blevel;
+
+// Face Barycenters (only needed for display)
+Eigen::MatrixXd B;
+
+// percentage of constrained faces
+double constraint_percentage = 0.002;
+
+// Random length factor
+double rand_factor = 5;
+
+// The set of parameters for calculating the curl-free fields
+igl::integrable_polyvector_fields_parameters params;
+
+// Solver data (needed for precomputation)
+igl::IntegrableFieldSolverData<Eigen::MatrixXd, Eigen::MatrixXi, Eigen::MatrixXd, Eigen::MatrixXd> ipfdata;
+
+//texture image
+Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic> texture_R, texture_G, texture_B;
+
+int display_mode = 1;
+
+int iter = 0;
+
+// Create a texture that hides the integer translation in the parametrization
+void line_texture(Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic> &texture_R,
+                  Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic> &texture_G,
+                  Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic> &texture_B)
+{
+  unsigned size = 128;
+  unsigned size2 = size/2;
+  unsigned lineWidth = 3;
+  texture_R.setConstant(size, size, 255);
+  for (unsigned i=0; i<size; ++i)
+    for (unsigned j=size2-lineWidth; j<=size2+lineWidth; ++j)
+      texture_R(i,j) = 0;
+  for (unsigned i=size2-lineWidth; i<=size2+lineWidth; ++i)
+    for (unsigned j=0; j<size; ++j)
+      texture_R(i,j) = 0;
+
+  texture_G = texture_R;
+  texture_B = texture_R;
+}
+
+// Create a random set of tangent vectors
+void generate_constraints()
+{
+  b.resize(42,1);
+  b<<
+  663,
+  513,
+  3872,
+  2601,
+  3549,
+  2721,
+  3796,
+  594,
+  868,
+  1730,
+  1581,
+  3081,
+  1471,
+  1650,
+  454,
+  2740,
+  2945,
+  3808,
+  3679,
+  3589,
+  450,
+  2656,
+  1791,
+  1792,
+  2917,
+  3744,
+  1536,
+  2809,
+  3866,
+  3658,
+  1665,
+  2670,
+  1601,
+  1793,
+  3614,
+  524,
+  2877,
+  449,
+  455,
+  3867,
+  3871,
+  2592;
+
+blevel.setOnes(b.rows(),1);
+bc.resize(b.rows(),6);
+bc<<
+-0.88046298335147721,0.27309862654264377,0.38755912468723769,-0.350632259447135,-0.92528970817792766,-0.14455440005410564,
+0.91471470003012889,0.392936119054695,-0.094330397492144599,0.32234487030777614,-0.85027369799342767,-0.41608703787410195,
+0.94335566040683105,0.073867667925654024,-0.32345581709658111,0.19950360079371404,-0.90525435056476755,0.37511714710727789,
+-0.92054671613540229,0.15077598183983737,0.36036141124232496,-0.27998315313687211,-0.89796618385425386,-0.33950871360506074,
+0.88944399663239937,0.23035525634795684,-0.39474780902172396,0.27297422303141039,-0.96047177712172194,0.054580572670497061,
+-0.83112706922096102,-0.55599928943162547,0.0096221078617792517,0.52546831822766438,-0.79091522174894457,-0.31358596675362826,
+0.90724658517664569,-0.41046292080872998,-0.091781394228251156,-0.34252813327252363,-0.84767620917196618,0.40511667741613094,
+-0.8932101465465786,0.23975524191487588,0.38038540729184012,-0.33645713296414853,-0.91759364410558497,-0.21170380718016926,
+-0.87839308390284521,0.27039404931731387,0.39409725734320344,-0.29712518405497651,-0.95481177255558192,-0.0071487054467167244,
+-0.91448048788760505,-0.17055891298176448,0.36692655188106316,0.29811257890714044,-0.89715315396744022,0.32595261714489804,
+0.82798126471567035,-0.56074230404745851,0.003885065171440813,-0.53510484459763941,-0.78801608401899037,0.30445600111594384,
+-0.87831929581593793,0.25312706437601257,0.40556368658667746,-0.26531767440854004,-0.9637845762158106,0.026941089342378936,
+-0.87482003689209031,-0.27011021313654948,0.4021571531272935,0.32303198334357713,-0.94388894366288889,0.0687313594225408,
+0.87408456883093666,-0.48487387939766946,-0.029554823793924323,-0.43846604347950752,-0.81368808449189478,0.38165328489525413,
+0.94988212941972827,-0.041936960956176939,-0.30978255521381903,-0.16637246118261648,-0.90677959514398765,-0.3873899456497869,
+0.87516493768857484,-0.27181042881473483,-0.40025669591913515,-0.36755520380602424,-0.91147911093961553,-0.18468622708756641,
+-0.87064244687577641,0.27922257819020818,0.40498948323008854,-0.32176729617260508,-0.94599403842079244,-0.039510585747255161,
+-0.91274615133859638,-0.1832657459904497,0.36511385835536858,0.29782083933521708,-0.91026141603074595,0.28762284704690655,
+0.875611546674125,0.28258715176515403,-0.39172556846369444,0.36000975242683186,-0.92250014843287764,0.13923524804764778,
+0.76763693171870195,-0.64088483679642994,0.00040868803559811201,-0.63058113310112196,-0.75518119878562417,0.17907761327874747,
+0.963064265211517,0.17044712473620016,-0.20845862597111031,0.061832174999749308,-0.89345471128813481,-0.44487690546019126,
+-0.88228998376691692,-0.46837234310148523,0.046815945597605227,0.41604986062280985,-0.82249303168905052,-0.38782434980116298,
+-0.96608602970701829,0.11121907649833783,0.23304098400879364,0.010641270548624404,-0.88457418950525291,0.46627810008860171,
+-0.96329451047686887,0.055809409647239343,0.26258140810033831,0.07182051046944142,-0.88891411988025926,0.45240855623364673,
+-0.71244584326772997,-0.70122065397026967,-0.026655484539588895,0.70046172163981768,-0.70836773631021255,-0.086997279682342638,
+0.88646445996853696,0.2549240118236365,-0.38625705094979518,0.35132981358631576,-0.91395520354543514,0.20310895597591658,
+-0.86109327343809683,-0.30822574449366841,0.40437020769461601,0.37896596246993836,-0.91928725525816557,0.10628142645421024,
+0.86443027504389158,-0.29669958642983363,-0.40586901212079213,-0.37200560813855077,-0.92052106924988175,-0.11938504337027039,
+0.95370728000967508,-0.24689991217686594,-0.17170572915195079,-0.14736898596800915,-0.88138264597997584,0.4488284898935197,
+-0.81439393313167019,0.57995723960933832,0.020300785774083896,-0.55494919604589421,-0.78855235001798585,0.26498411478639117,
+0.89527216270596455,0.22395367264061938,-0.38513959442592183,0.33680943342191538,-0.90609008785063272,0.25604717974787594,
+-0.9003647006267198,0.20802946062196581,0.38218732236782926,-0.32431023000528064,-0.90640636884236436,-0.27064805418831556,
+-0.87050937437709508,-0.28614105672408718,0.40042068475344922,0.37746788793940733,-0.91025870352880611,0.17013843253251126,
+-0.95715079751532439,0.0030851788865496879,0.28957353554324744,0.12908381923211401,-0.89056292562302997,0.43615942397041058,
+-0.87324677619075319,-0.28591869514051466,0.39457644080913162,0.3438918663433696,-0.93530416305293196,0.083333707698197687,
+0.91999856277124803,-0.1621255206103257,-0.35681642348085474,-0.27672206872177485,-0.91342693749618353,-0.2984562389005877,
+-0.8669467282645521,0.29036174243712859,0.40508447128995645,-0.34873789125620602,-0.93406639205959408,-0.07682355385522964,
+-0.9570365266718136,-0.22821899053183647,0.17887755302603078,0.12590409644663494,-0.88275887883510706,-0.45264215483728532,
+-0.94033215083998489,0.087395510869996196,0.32884262311388451,-0.2131320783418921,-0.90465024471116184,-0.36902933748646671,
+-0.96131014054749453,0.18866284908038999,0.20072155603578434,-0.08260532909072589,-0.89255302833360861,0.44331191188407898,
+-0.95240414686152941,-0.02752900142620229,0.30359264668538755,0.15128346580527452,-0.9073021943457209,0.39232134929083828,
+-0.94070423353276911,-0.31552769387286655,0.12457053990729766,0.22959741970407915,-0.86253407908715607,-0.45091017650802745;
+
+}
+
+void drawCuts(igl::viewer::Viewer& viewer,
+              const Eigen::MatrixXi &cuts)
+{
+  int maxCutNum = cuts.sum();
+  Eigen::MatrixXd start(maxCutNum,3);
+  Eigen::MatrixXd end(maxCutNum,3);
+  int ind = 0;
+  for (unsigned int i=0;i<F.rows();i++)
+    for (int j=0;j<3;j++)
+      if (cuts(i,j))
+      {
+        start.row(ind) = V.row(F(i,j));
+        end.row(ind) = V.row(F(i,(j+1)%3));
+        ind++;
+      }
+  viewer.data.add_edges(start, end , Eigen::RowVector3d(1.,0,1.));
+}
+
+void drawField(igl::viewer::Viewer &viewer,
+               const Eigen::MatrixXd &field)
+{
+  for (int n=0; n<2; ++n)
+  {
+    Eigen::MatrixXd VF = field.block(0,n*3,F.rows(),3);
+    Eigen::VectorXd c = VF.rowwise().norm();
+    Eigen::MatrixXd C2;
+    igl::jet(c,1,1+rand_factor,C2);
+    viewer.data.add_edges(B - global_scale*VF, B + global_scale*VF , C2);
+  }
+}
+
+void drawConstraints(igl::viewer::Viewer &viewer)
+{
+  for (int n=0; n<2; ++n)
+  {
+    Eigen::MatrixXd Bc = igl::slice(B, b, 1);
+    Eigen::MatrixXd color;
+    color.setZero(b.rows(),3);
+    color.col(2).setOnes();
+    for (int i =0; i<b.rows(); ++i)
+      if (blevel[i] ==1 && n>0)
+        color.row(i)<<0.7,0.7,0.7;
+    // Eigen::RowVector3d color; color<<0.5,0.5,0.5;
+    viewer.data.add_edges(Bc - global_scale*bc.block(0,n*3,bc.rows(),3), Bc + global_scale*bc.block(0,n*3,bc.rows(),3) , color);
+  }
+
+}
+
+
+void colorEdgeMeshFaces(const Eigen::VectorXd &values,
+                        const double &minimum,
+                        const double &maximum,
+                        Eigen::MatrixXd &C)
+{
+  C.setConstant(Fbs.rows(),3,1);
+
+  Eigen::MatrixXd colors;
+  igl::jet(values, minimum, maximum, colors);
+
+  for (int ei = 0; ei<E.rows(); ++ei)
+  {
+    const Eigen::RowVector3d &this_color = colors.row(ei);
+    int f0 = E2F(ei,0);
+    int f1 = E2F(ei,1);
+    if(f0 != -1)
+    {
+      int i0 = -1;
+      for (int k = 0; k<3; k++)
+        if (F2E(f0,k)== ei)
+        {
+          i0 = k;
+          break;
+        }
+      C.row(3*f0+i0) = this_color;
+    }
+    if(f1 != -1)
+    {
+      int i1 = -1;
+      for (int k = 0; k<3; k++)
+        if (F2E(f1,k)== ei)
+        {
+          i1 = k;
+          break;
+        }
+      C.row(3*f1+i1) = this_color;
+    }
+  }
+
+}
+
+void update_display(igl::viewer::Viewer& viewer)
+{
+  using namespace std;
+  using namespace Eigen;
+
+  viewer.data.clear();
+  viewer.data.lines.resize(0,9);
+  viewer.data.points.resize(0,6);
+  viewer.core.show_texture = false;
+
+  if (display_mode == 1)
+  {
+    cerr<< "Displaying original field, its singularities and its cuts"  <<endl;
+
+    viewer.data.set_mesh(V, F);
+
+    // Highlight in red the constrained faces
+    MatrixXd C = MatrixXd::Constant(F.rows(),3,1);
+    for (unsigned i=0; i<b.size();++i)
+      C.row(b(i)) << 1, 0, 0;
+    viewer.data.set_colors(C);
+
+    //Draw constraints
+    drawConstraints(viewer);
+
+    // Draw Field
+    drawField(viewer,two_pv_ori);
+
+    // Draw Cuts
+    drawCuts(viewer,cuts_ori);
+
+    //Draw Singularities
+    Eigen::MatrixXd singular_points = igl::slice(V, singularities_ori, 1);
+    viewer.data.add_points(singular_points,Eigen::RowVector3d(239./255.,205./255.,57./255.));
+
+  }
+
+  if (display_mode == 2)
+  {
+    cerr<< "Displaying current field, its singularities and its cuts"  <<endl;
+
+    viewer.data.set_mesh(V, F);
+
+    // Highlight in red the constrained faces
+    MatrixXd C = MatrixXd::Constant(F.rows(),3,1);
+    for (unsigned i=0; i<b.size();++i)
+      C.row(b(i)) << 1, 0, 0;
+    viewer.data.set_colors(C);
+
+    //Draw constraints
+    drawConstraints(viewer);
+
+    // Draw Field
+    drawField(viewer,two_pv);
+
+    // Draw Cuts
+    drawCuts(viewer,cuts);
+
+    //Draw Singularities
+    Eigen::MatrixXd singular_points = igl::slice(V, singularities, 1);
+    viewer.data.add_points(singular_points,Eigen::RowVector3d(239./255.,205./255.,57./255.));
+  }
+
+  if (display_mode == 3)
+  {
+    cerr<< "Displaying original field and its curl"  <<endl;
+
+    viewer.data.set_mesh(Vbs, Fbs);
+    Eigen::MatrixXd C;
+    colorEdgeMeshFaces(curl_ori, 0, 0.2, C);
+    viewer.data.set_colors(C);
+
+    // Draw Field
+    drawField(viewer,two_pv_ori);
+
+  }
+
+  if (display_mode == 4)
+  {
+    cerr<< "Displaying current field and its curl"  <<endl;
+
+    viewer.data.set_mesh(Vbs, Fbs);
+    Eigen::MatrixXd C;
+    colorEdgeMeshFaces(curl, 0, 0.2, C);
+    viewer.data.set_colors(C);
+
+    // Draw Field
+    drawField(viewer,two_pv);
+  }
+
+  if (display_mode == 5)
+  {
+    cerr<< "Displaying original poisson-integrated field and original poisson error"  <<endl;
+
+    viewer.data.set_mesh(V, F);
+    Eigen::MatrixXd C;
+    igl::jet(poisson_error_ori, 0, 0.5, C);
+    viewer.data.set_colors(C);
+
+    // Draw Field
+    drawField(viewer,two_pv_poisson_ori);
+  }
+
+  if (display_mode == 6)
+  {
+    cerr<< "Displaying current poisson-integrated field and current poisson error"  <<endl;
+
+    viewer.data.set_mesh(V, F);
+    Eigen::MatrixXd C;
+    igl::jet(poisson_error, 0, 0.5, C);
+    viewer.data.set_colors(C);
+
+    // Draw Field
+    drawField(viewer,two_pv_poisson);
+  }
+
+  if (display_mode == 7)
+  {
+    cerr<< "Displaying original texture with cuts and singularities"  <<endl;
+
+    viewer.data.set_mesh(V, F);
+    MatrixXd C = MatrixXd::Constant(F.rows(),3,1);
+    viewer.data.set_colors(C);
+    viewer.data.set_uv(uv_scale*scalars_ori, Fcut_ori);
+    viewer.data.set_texture(texture_R, texture_B, texture_G);
+    viewer.core.show_texture = true;
+
+    // Draw Cuts
+    drawCuts(viewer,cuts_ori);
+
+    //Draw Singularities
+    Eigen::MatrixXd singular_points = igl::slice(V, singularities_ori, 1);
+    viewer.data.add_points(singular_points,Eigen::RowVector3d(239./255.,205./255.,57./255.));
+
+  }
+  if (display_mode == 8)
+  {
+    cerr<< "Displaying current texture with cuts and singularities"  <<endl;
+
+    viewer.data.set_mesh(V, F);
+    MatrixXd C = MatrixXd::Constant(F.rows(),3,1);
+    viewer.data.set_colors(C);
+    viewer.data.set_uv(uv_scale*scalars, Fcut);
+    viewer.data.set_texture(texture_R, texture_B, texture_G);
+    viewer.core.show_texture = true;
+
+    // Draw Cuts
+    drawCuts(viewer,cuts);
+
+    //Draw Singularities
+    Eigen::MatrixXd singular_points = igl::slice(V, singularities, 1);
+    viewer.data.add_points(singular_points,Eigen::RowVector3d(239./255.,205./255.,57./255.));
+
+  }
+
+  if (display_mode == 9)
+  {
+    cerr<< "Displaying original field overlayed onto the current integrated field"  <<endl;
+
+    viewer.data.set_mesh(V, F);
+
+    // Highlight in red the constrained faces
+    MatrixXd C = MatrixXd::Constant(F.rows(),3,1);
+    for (unsigned i=0; i<b.size();++i)
+      C.row(b(i)) << 1, 0, 0;
+    viewer.data.set_colors(C);
+
+    // Draw Field
+    drawField(viewer,two_pv_ori);
+
+    // Draw Integrated Field
+    drawField(viewer,two_pv_poisson_ori);
+
+  }
+
+  if (display_mode == 0)
+  {
+    cerr<< "Displaying current field overlayed onto the current integrated field"  <<endl;
+
+    viewer.data.set_mesh(V, F);
+
+    // Highlight in red the constrained faces
+    MatrixXd C = MatrixXd::Constant(F.rows(),3,1);
+    for (unsigned i=0; i<b.size();++i)
+      C.row(b(i)) << 1, 0, 0;
+    viewer.data.set_colors(C);
+
+    // Draw Field
+    drawField(viewer,two_pv);
+
+    // Draw Integrated Field
+    drawField(viewer,two_pv_poisson);
+  }
+
+}
+
+bool key_down(igl::viewer::Viewer& viewer, unsigned char key, int modifier)
+{
+
+  if (key == '1')
+  {
+    display_mode = 1;
+    update_display(viewer);
+  }
+
+  if (key == '2')
+  {
+    display_mode = 2;
+    update_display(viewer);
+  }
+
+  if (key == '3')
+  {
+    display_mode = 3;
+    update_display(viewer);
+  }
+
+  if (key == '4')
+  {
+    display_mode = 4;
+    update_display(viewer);
+  }
+
+  if (key == '5')
+  {
+    display_mode = 5;
+    update_display(viewer);
+  }
+
+  if (key == '6')
+  {
+    display_mode = 6;
+    update_display(viewer);
+  }
+
+  if (key == '7')
+  {
+    display_mode = 7;
+    update_display(viewer);
+  }
+
+  if (key == '8')
+  {
+    display_mode = 8;
+    update_display(viewer);
+  }
+
+  if (key == '9')
+  {
+    display_mode = 9;
+    update_display(viewer);
+  }
+
+  if (key == '0')
+  {
+    display_mode = 0;
+    update_display(viewer);
+  }
+
+  if (key == 'A')
+  {
+    //do a batch of iterations
+    printf("--Improving Curl--\n");
+    for (int bi = 0; bi<5; ++bi)
+    {
+      printf("\n\n **** Batch %d ****\n", iter);
+      igl::integrable_polyvector_fields_solve(ipfdata, params, two_pv, iter ==0);
+      iter++;
+      params.wSmooth *= params.redFactor_wsmooth;
+    }
+    // Post process current field
+    // Compute curl_minimizing matchings and curl
+    printf("--Matchings and curl--\n");
+    Eigen::MatrixXi match_ab, match_ba;  // matchings across interior edges
+    double avgCurl = igl::polyvector_field_matchings(two_pv, V, F, true, match_ab, match_ba, curl);
+    double maxCurl = curl.maxCoeff();
+    printf("curl -- max: %.5g, avg: %.5g\n", maxCurl,  avgCurl);
+    // Compute singularities
+    printf("--Singularities--\n");
+    igl::polyvector_field_singularities_from_matchings(V, F, match_ab, match_ba, singularities);
+    printf("#singularities: %ld\n", singularities.rows());
+    // Get mesh cuts based on singularities
+    printf("--Cuts--\n");
+    igl::polyvector_field_cut_mesh_with_singularities(V, F, singularities, cuts);
+    // Comb field
+    printf("--Combing--\n");
+    Eigen::MatrixXd combed;
+    igl::polyvector_field_comb_from_matchings_and_cuts(V, F, two_pv, match_ab, match_ba, cuts, combed);
+    // Reconstruct integrable vector fields from combed field
+    printf("--Cut mesh--\n");
+    igl::cut_mesh(V, F, cuts, Vcut, Fcut);
+    printf("--Poisson--\n");
+    double avgPoisson = igl::polyvector_field_poisson_reconstruction(Vcut, Fcut, combed, scalars, two_pv_poisson, poisson_error);
+    double maxPoisson = poisson_error.maxCoeff();
+    printf("poisson error -- max: %.5g, avg: %.5g\n", maxPoisson, avgPoisson);
+
+    update_display(viewer);
+  }
+
+  return false;
+}
+
+int main(int argc, char *argv[])
+{
+
+  // Load a mesh
+  igl::readOBJ("../shared/inspired_mesh.obj", V, F);
+
+  printf("--Initialization--\n");
+  V_border = igl::is_border_vertex(V,F);
+  igl::adjacency_list(F, VV);
+  igl::vertex_triangle_adjacency(V,F,VF,VFi);
+  igl::triangle_triangle_adjacency(V,F,TT,TTi);
+  igl::edge_topology(V,F,E,F2E,E2F);
+
+  // Generate "subdivided" mesh for visualization of curl terms
+  igl::false_barycentric_subdivision(V, F, Vbs, Fbs);
+
+  // Compute scale for visualizing fields
+  global_scale =  .2*igl::avg_edge_length(V, F);
+
+  //Compute scale for visualizing texture
+  uv_scale = 0.6/igl::avg_edge_length(V, F);
+
+  // Compute face barycenters
+  igl::barycenter(V, F, B);
+
+  // Compute local basis for faces
+  igl::local_basis(V,F,B1,B2,B3);
+
+  //Generate random vectors for constraints
+  generate_constraints();
+
+  // Interpolate a 2-PolyVector field to be used as the original field
+  printf("--Initial solution--\n");
+  igl::n_polyvector(V, F, b, bc, two_pv_ori);
+
+  // Post process original field
+  // Compute curl_minimizing matchings and curl
+  Eigen::MatrixXi match_ab, match_ba;  // matchings across interior edges
+  printf("--Matchings and curl--\n");
+  double avgCurl = igl::polyvector_field_matchings(two_pv_ori, V, F, true, match_ab, match_ba, curl_ori);
+  double maxCurl = curl_ori.maxCoeff();
+  printf("original curl -- max: %.5g, avg: %.5g\n", maxCurl,  avgCurl);
+
+  printf("--Singularities--\n");
+  // Compute singularities
+  igl::polyvector_field_singularities_from_matchings(V, F, V_border, VF, TT, E2F, F2E, match_ab, match_ba, singularities_ori);
+  printf("original #singularities: %ld\n", singularities.rows());
+
+  printf("--Cuts--\n");
+ // Get mesh cuts based on singularities
+  igl::polyvector_field_cut_mesh_with_singularities(V, F, VF, VV, TT, TTi, singularities_ori, cuts_ori);
+
+  printf("--Combing--\n");
+// Comb field
+  Eigen::MatrixXd combed;
+  igl::polyvector_field_comb_from_matchings_and_cuts(V, F, TT, E2F, F2E, two_pv_ori, match_ab, match_ba, cuts_ori, combed);
+
+  printf("--Cut mesh--\n");
+  // Reconstruct integrable vector fields from combed field
+  igl::cut_mesh(V, F, VF, VFi, TT, TTi, V_border, cuts_ori, Vcut_ori, Fcut_ori);
+
+  printf("--Poisson--\n");
+  double avgPoisson = igl::polyvector_field_poisson_reconstruction(Vcut_ori, Fcut_ori, combed, scalars_ori, two_pv_poisson_ori, poisson_error_ori);
+  double maxPoisson = poisson_error_ori.maxCoeff();
+  printf("poisson error -- max: %.5g, avg: %.5g\n", maxPoisson, avgPoisson);
+
+
+  // Set the curl-free 2-PolyVector to equal the original field
+  two_pv = two_pv_ori;
+  singularities = singularities_ori;
+  curl = curl_ori;
+  cuts = cuts_ori;
+  two_pv_poisson = two_pv_poisson_ori;
+  poisson_error = poisson_error_ori;
+  Vcut = Vcut_ori;
+  Fcut = Fcut_ori;
+  scalars = scalars_ori;
+
+  printf("--Integrable - Precomputation--\n");
+  // Precompute stuff for solver
+  igl::integrable_polyvector_fields_precompute(V, F, b, bc, blevel, two_pv_ori, ipfdata);
+
+  cerr<<"Done. Press keys 1-0 for various visualizations, 'A' to improve integrability." <<endl;
+
+  igl::viewer::Viewer viewer;
+  viewer.callback_key_down = &key_down;
+  viewer.core.show_lines = false;
+  key_down(viewer,'2',0);
+
+  // Replace the standard texture with an integer shift invariant texture
+  line_texture(texture_R, texture_G, texture_B);
+
+  viewer.launch();
+
+  return 0;
+}
diff --git a/tutorial/CMakeLists.txt b/tutorial/CMakeLists.txt
index 612da5873..fe257438a 100644
--- a/tutorial/CMakeLists.txt
+++ b/tutorial/CMakeLists.txt
@@ -122,3 +122,4 @@ add_subdirectory("701_Statistics")
 add_subdirectory("702_WindingNumber")
 add_subdirectory("703_Decimation")
 add_subdirectory("704_SignedDistance")
+add_subdirectory("705_Integrable")