compilation fixes for VS 2015

Former-commit-id: c50570af1c1954b270d70c1e70bb82f05ed7c52a

include/igl/angle_bound_frame_fields.cpp

@@ -38,7 +38,7 @@ namespace igl {
       Eigen::VectorXi indInteriorToFull;
       Eigen::VectorXi indFullToInterior;
 
-#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
+//#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
       Eigen::PlainObjectBase<DerivedV> B1, B2, FN;
 
       //laplacians

include/igl/colon.cpp

@@ -1,9 +1,9 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-// 
+//
 // Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
-// 
-// This Source Code Form is subject to the terms of the Mozilla Public License 
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "colon.h"
 
@@ -11,9 +11,9 @@
 
 template <typename L,typename S,typename H,typename T>
 IGL_INLINE void igl::colon(
-  const L low, 
-  const S step, 
-  const H hi, 
+  const L low,
+  const S step,
+  const H hi,
   Eigen::Matrix<T,Eigen::Dynamic,1> & I)
 {
   const int size = ((hi-low)/step)+1;
@@ -22,8 +22,8 @@ IGL_INLINE void igl::colon(
 
 template <typename L,typename H,typename T>
 IGL_INLINE void igl::colon(
-  const L low, 
-  const H hi, 
+  const L low,
+  const H hi,
   Eigen::Matrix<T,Eigen::Dynamic,1> & I)
 {
   return igl::colon(low,(T)1,hi,I);
@@ -31,7 +31,7 @@ IGL_INLINE void igl::colon(
 
 template <typename T,typename L,typename H>
 IGL_INLINE Eigen::Matrix<T,Eigen::Dynamic,1> igl::colon(
-  const L low, 
+  const L low,
   const H hi)
 {
   Eigen::Matrix<T,Eigen::Dynamic,1> I;
@@ -53,4 +53,5 @@ template void igl::colon<int, int, int>(int, int, Eigen::Matrix<int, -1, 1, 0, -
 template void igl::colon<int,long long int,int>(int,long long int,Eigen::Matrix<int,-1,1,0,-1,1> &);
 template void igl::colon<int, int, int, int>(int, int, int, Eigen::Matrix<int, -1, 1, 0, -1, 1>&);
 template void igl::colon<int, long, long>(int, long, Eigen::Matrix<long, -1, 1, 0, -1, 1>&);
+template void igl::colon<int, long long,long>(int, long long, class Eigen::Matrix<long,-1,1,0,-1,1> &);
 #endif

include/igl/comb_cross_field.cpp

@@ -27,7 +27,7 @@ namespace igl {
     const Eigen::PlainObjectBase<DerivedF> &F;
     const Eigen::PlainObjectBase<DerivedV> &PD1;
     const Eigen::PlainObjectBase<DerivedV> &PD2;
-#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
+//#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
     Eigen::PlainObjectBase<DerivedV> N;
 
   private:

include/igl/comb_line_field.cpp

@@ -25,12 +25,12 @@ public:
     const Eigen::PlainObjectBase<DerivedV> &V;
     const Eigen::PlainObjectBase<DerivedF> &F;
     const Eigen::PlainObjectBase<DerivedV> &PD1;
-#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
+//#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
     Eigen::PlainObjectBase<DerivedV> N;
 
 private:
     // internal
-#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
+//#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
     Eigen::PlainObjectBase<DerivedF> TT;
     Eigen::PlainObjectBase<DerivedF> TTi;
 

include/igl/copyleft/quadprog.cpp

@@ -6,7 +6,7 @@
 // v. 2.0. If a copy of the MPL was not distributed with this file, You can 
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "quadprog.h"
-
+#include <vector>
 /*
  FILE eiquadprog.hh
  
@@ -307,7 +307,7 @@ IGL_INLINE bool igl::copyleft::quadprog(
   VectorXi A(m + p), A_old(m + p), iai(m + p);
   int q;
   int iq, iter = 0;
-  bool iaexcl[m + p];
+  std::vector<bool> iaexcl(m + p);
  	
   me = p; /* number of equality constraints */
   mi = m; /* number of inequality constraints */

include/igl/cross_field_missmatch.cpp

@@ -28,7 +28,7 @@ namespace igl {
     const Eigen::PlainObjectBase<DerivedF> &F;
     const Eigen::PlainObjectBase<DerivedV> &PD1;
     const Eigen::PlainObjectBase<DerivedV> &PD2;
-#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
+//#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
     Eigen::PlainObjectBase<DerivedV> N;
 
   private:
@@ -36,7 +36,7 @@ namespace igl {
     std::vector<bool> V_border; // bool
     std::vector<std::vector<int> > VF;
     std::vector<std::vector<int> > VFi;
-#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
+//#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
     Eigen::PlainObjectBase<DerivedF> TT;
     Eigen::PlainObjectBase<DerivedF> TTi;
 

include/igl/cut_mesh.cpp

@@ -29,7 +29,7 @@ namespace igl {
     int num_scalar_variables;
 
     // per face indexes of vertex in the solver
-#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
+//#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
     Eigen::PlainObjectBase<DerivedF> HandleS_Index;
 
     // per vertex variable indexes

include/igl/cut_mesh_from_singularities.cpp

@@ -30,7 +30,7 @@ namespace igl {
     const Eigen::PlainObjectBase<DerivedM> &Handle_MMatch;
 
     Eigen::VectorXi F_visited;
-#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
+//#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
     Eigen::PlainObjectBase<DerivedF> TT;
     Eigen::PlainObjectBase<DerivedF> TTi;
 

include/igl/eigs.cpp

@@ -153,4 +153,5 @@ IGL_INLINE bool igl::eigs(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 template bool igl::eigs<double, double, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::SparseMatrix<double, 0, int> const&, Eigen::SparseMatrix<double, 0, int> const&, unsigned long, igl::EigsType, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template bool igl::eigs<double, double, Eigen::Matrix<double,-1,-1,0,-1,-1>, Eigen::Matrix<double,-1,1,0,-1,1> >(Eigen::SparseMatrix<double,0,int> const &,Eigen::SparseMatrix<double,0,int> const &,unsigned long long, igl::EigsType, Eigen::PlainObjectBase< Eigen::Matrix<double,-1,-1,0,-1,-1> > &, Eigen::PlainObjectBase<Eigen::Matrix<double,-1,1,0,-1,1> > &);
 #endif

include/igl/line_field_missmatch.cpp

@@ -31,7 +31,7 @@ public:
     const Eigen::PlainObjectBase<DerivedF> &F;
     const Eigen::PlainObjectBase<DerivedV> &PD1;
     const Eigen::PlainObjectBase<DerivedV> &PD2;
-#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
+//#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
     Eigen::PlainObjectBase<DerivedV> N;
 
 private:

include/igl/list_to_matrix.cpp

@@ -135,4 +135,5 @@ template bool igl::list_to_matrix<bool, Eigen::Array<bool, -1, -1, 0, -1, -1> >(
 template bool igl::list_to_matrix<bool, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(std::vector<bool, std::allocator<bool> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template bool igl::list_to_matrix<bool, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(std::vector<std::vector<bool, std::allocator<bool> >, std::allocator<std::vector<bool, std::allocator<bool> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template bool igl::list_to_matrix<double, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(std::vector<double, std::allocator<double> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template bool igl::list_to_matrix<unsigned long long,class Eigen::Matrix<int,-1,1,0,-1,1> >(class std::vector<unsigned long long,class std::allocator<unsigned long long> > const &,class Eigen::PlainObjectBase<class Eigen::Matrix<int,-1,1,0,-1,1> > &);
 #endif

include/igl/n_polyvector.cpp

@@ -39,7 +39,7 @@ namespace igl {
     Eigen::VectorXi indInteriorToFull;
     Eigen::VectorXi indFullToInterior;
 
-#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
+//#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
     Eigen::PlainObjectBase<DerivedV> B1, B2, FN;
 
     IGL_INLINE void computek();

include/igl/n_polyvector_general.cpp

@@ -7,468 +7,465 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 
 #include <igl/n_polyvector_general.h>
-#include <igl/edge_topology.h>
-#include <igl/local_basis.h>
-#include <igl/nchoosek.h>
-#include <igl/slice.h>
-#include <igl/polyroots.h>
-#include <Eigen/Sparse>
-#include <Eigen/Geometry>
-#include <iostream>
-
-namespace igl {
-  template <typename DerivedV, typename DerivedF>
-  class GeneralPolyVectorFieldFinder
-  {
-  private:
-    const Eigen::PlainObjectBase<DerivedV> &V;
-    const Eigen::PlainObjectBase<DerivedF> &F; int numF;
-    const int n;
-
-    Eigen::MatrixXi EV; int numE;
-    Eigen::MatrixXi F2E;
-    Eigen::MatrixXi E2F;
-    Eigen::VectorXd K;
-
-    Eigen::VectorXi isBorderEdge;
-    int numInteriorEdges;
-    Eigen::Matrix<int,Eigen::Dynamic,2> E2F_int;
-    Eigen::VectorXi indInteriorToFull;
-    Eigen::VectorXi indFullToInterior;
-
-#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
-    Eigen::PlainObjectBase<DerivedV> B1, B2, FN;
-
-    IGL_INLINE void computek();
-    IGL_INLINE void setFieldFromGeneralCoefficients(const  std::vector<Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1>> &coeffs,
-                                                    std::vector<Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 2> > &pv);
-    IGL_INLINE void computeCoefficientLaplacian(int n, Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &D);
-    IGL_INLINE void getGeneralCoeffConstraints(const Eigen::VectorXi &isConstrained,
-                                    const Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> &cfW,
-                                    int k,
-                                    const Eigen::VectorXi &rootsIndex,
-                                    Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1> &Ck);
-    IGL_INLINE void precomputeInteriorEdges();
-
-
-    IGL_INLINE void minQuadWithKnownMini(const Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &Q,
-                                         const Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &f,
-                                         const Eigen::VectorXi isConstrained,
-                                         const Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> &xknown,
-                                         Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> &x);
-
-  public:
-    IGL_INLINE GeneralPolyVectorFieldFinder(const Eigen::PlainObjectBase<DerivedV> &_V,
-                                     const Eigen::PlainObjectBase<DerivedF> &_F,
-                                     const int &_n);
-    IGL_INLINE bool solve(const Eigen::VectorXi &isConstrained,
-               const Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> &cfW,
-               const Eigen::VectorXi &rootsIndex,
-               Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> &output);
-
-  };
-}
-
-template<typename DerivedV, typename DerivedF>
-IGL_INLINE igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::
-          GeneralPolyVectorFieldFinder(const Eigen::PlainObjectBase<DerivedV> &_V,
-                                const Eigen::PlainObjectBase<DerivedF> &_F,
-                                const int &_n):
-V(_V),
-F(_F),
-numF(_F.rows()),
-n(_n)
-{
-
-  igl::edge_topology(V,F,EV,F2E,E2F);
-  numE = EV.rows();
-
-
-  precomputeInteriorEdges();
-
-  igl::local_basis(V,F,B1,B2,FN);
-
-  computek();
-
-};
-
-
-template<typename DerivedV, typename DerivedF>
-IGL_INLINE void igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::
-precomputeInteriorEdges()
-{
-  // Flag border edges
-  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)))
-      isBorderEdge[i] = 1;
-      else
-      {
-        indFullToInterior[i] = numInteriorEdges;
-        numInteriorEdges++;
-      }
-  }
-
-  E2F_int.resize(numInteriorEdges, 2);
-  indInteriorToFull.setZero(numInteriorEdges,1);
-  int ii = 0;
-  for (int k=0; k<numE; ++k)
-  {
-    if (isBorderEdge[k])
-      continue;
-    E2F_int.row(ii) = E2F.row(k);
-    indInteriorToFull[ii] = k;
-    ii++;
-  }
-
-}
-
-
-template<typename DerivedV, typename DerivedF>
-IGL_INLINE void igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::
-minQuadWithKnownMini(const Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &Q,
-                          const Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &f,
-                     const Eigen::VectorXi isConstrained,
-                          const Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> &xknown,
-                          Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> &x)
-{
-  int N = Q.rows();
-
-  int nc = xknown.rows();
-  Eigen::VectorXi known; known.setZero(nc,1);
-  Eigen::VectorXi unknown; unknown.setZero(N-nc,1);
-
-  int indk = 0, indu = 0;
-  for (int i = 0; i<N; ++i)
-    if (isConstrained[i])
-    {
-      known[indk] = i;
-      indk++;
-    }
-    else
-    {
-      unknown[indu] = i;
-      indu++;
-    }
-
-  Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar>> Quu, Quk;
-
-  igl::slice(Q,unknown, unknown, Quu);
-  igl::slice(Q,unknown, known, Quk);
-
-
-  std::vector<typename Eigen::Triplet<std::complex<typename DerivedV::Scalar> > > tripletList;
-
-  Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > fu(N-nc,1);
-
-  igl::slice(f,unknown, Eigen::VectorXi::Zero(1,1), fu);
-
-  Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > rhs = (Quk*xknown).sparseView()+.5*fu;
-
-  Eigen::SparseLU< Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar>>> solver;
-  solver.compute(-Quu);
-  if(solver.info()!=Eigen::Success)
-  {
-    std::cerr<<"Decomposition failed!"<<std::endl;
-    return;
-  }
-  Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar>>  b  = solver.solve(rhs);
-  if(solver.info()!=Eigen::Success)
-  {
-    std::cerr<<"Solving failed!"<<std::endl;
-    return;
-  }
-
-  indk = 0, indu = 0;
-  x.setZero(N,1);
-  for (int i = 0; i<N; ++i)
-    if (isConstrained[i])
-      x[i] = xknown[indk++];
-    else
-      x[i] = b.coeff(indu++,0);
-
-}
-
-
-
-template<typename DerivedV, typename DerivedF>
-IGL_INLINE bool igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::
-                     solve(const Eigen::VectorXi &isConstrained,
-                           const Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> &cfW,
-                           const Eigen::VectorXi &rootsIndex,
-                           Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> &output)
-{
-
-  // polynomial is of the form:
-  // z^(2n) +
-  // -c[0]z^(2n-1) +
-  // c[1]z^(2n-2) +
-  // -c[2]z^(2n-3) +
-  // ... +
-  // (-1)^n c[n-1]
-
-  std::vector<Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1>> coeffs(n,Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1>::Zero(numF, 1));
-
-  for (int i =0; i<n; ++i)
-  {
-    int degree = i+1;
-
-    Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1> Ck;
-    getGeneralCoeffConstraints(isConstrained,
-                               cfW,
-                               i,
-                               rootsIndex,
-                               Ck);
-
-    Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > DD;
-    computeCoefficientLaplacian(degree, DD);
-    Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > f; f.resize(numF,1);
-
-    if (isConstrained.sum() == numF)
-      coeffs[i] = Ck;
-    else
-      minQuadWithKnownMini(DD, f, isConstrained, Ck, coeffs[i]);
-  }
-
-  std::vector<Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 2> > pv;
-  setFieldFromGeneralCoefficients(coeffs, pv);
-
-  output.setZero(numF,3*n);
-  for (int fi=0; fi<numF; ++fi)
-  {
-    const Eigen::Matrix<typename DerivedV::Scalar, 1, 3> &b1 = B1.row(fi);
-    const Eigen::Matrix<typename DerivedV::Scalar, 1, 3> &b2 = B2.row(fi);
-    for (int i=0; i<n; ++i)
-      output.block(fi,3*i, 1, 3) = pv[i](fi,0)*b1 + pv[i](fi,1)*b2;
-  }
-  return true;
-}
-
-template<typename DerivedV, typename DerivedF>
-IGL_INLINE void igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::setFieldFromGeneralCoefficients(const  std::vector<Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1>> &coeffs,
-                                                            std::vector<Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 2>> &pv)
-{
-  pv.assign(n, Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 2>::Zero(numF, 2));
-  for (int i = 0; i <numF; ++i)
-  {
-
-    //    poly coefficients: 1, 0, -Acoeff, 0, Bcoeff
-    //    matlab code from roots (given there are no trailing zeros in the polynomial coefficients)
-    Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1> polyCoeff;
-    polyCoeff.setZero(n+1,1);
-    polyCoeff[0] = 1.;
-    int sign = 1;
-    for (int k =0; k<n; ++k)
-    {
-      sign = -sign;
-      int degree = k+1;
-      polyCoeff[degree] = (1.*sign)*coeffs[k](i);
-    }
-
-    Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1> roots;
-    igl::polyRoots<std::complex<typename DerivedV::Scalar>, typename DerivedV::Scalar >(polyCoeff,roots);
-    for (int k=0; k<n; ++k)
-    {
-      pv[k](i,0) = real(roots[k]);
-      pv[k](i,1) = imag(roots[k]);
-    }
-  }
-
-}
-
-
-template<typename DerivedV, typename DerivedF>
-IGL_INLINE void igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::computeCoefficientLaplacian(int n, Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &D)
-{
-  std::vector<Eigen::Triplet<std::complex<typename DerivedV::Scalar> >> tripletList;
-
-  // For every non-border edge
-  for (unsigned eid=0; eid<numE; ++eid)
-  {
-    if (!isBorderEdge[eid])
-    {
-      int fid0 = E2F(eid,0);
-      int fid1 = E2F(eid,1);
-
-      tripletList.push_back(Eigen::Triplet<std::complex<typename DerivedV::Scalar> >(fid0,
-                                           fid0,
-                                           std::complex<typename DerivedV::Scalar>(1.)));
-      tripletList.push_back(Eigen::Triplet<std::complex<typename DerivedV::Scalar> >(fid1,
-                                           fid1,
-                                           std::complex<typename DerivedV::Scalar>(1.)));
-      tripletList.push_back(Eigen::Triplet<std::complex<typename DerivedV::Scalar> >(fid0,
-                                           fid1,
-                                                                                     -1.*std::polar(1.,-1.*n*K[eid])));
-      tripletList.push_back(Eigen::Triplet<std::complex<typename DerivedV::Scalar> >(fid1,
-                                           fid0,
-                                                                                     -1.*std::polar(1.,1.*n*K[eid])));
-
-    }
-  }
-  D.resize(numF,numF);
-  D.setFromTriplets(tripletList.begin(), tripletList.end());
-
-
-}
-
-//this gives the coefficients without the (-1)^k that multiplies them
-template<typename DerivedV, typename DerivedF>
-IGL_INLINE void igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::getGeneralCoeffConstraints(const Eigen::VectorXi &isConstrained,
-                                                       const Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> &cfW,
-                                                       int k,
-                                                       const Eigen::VectorXi &rootsIndex,
-                                                       Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1> &Ck)
-{
-  int numConstrained = isConstrained.sum();
-  Ck.resize(numConstrained,1);
-  // int n = rootsIndex.cols();
-
-  Eigen::MatrixXi allCombs;
-  {
-    Eigen::VectorXi V = Eigen::VectorXi::LinSpaced(n,0,n-1);
-    igl::nchoosek(V,k+1,allCombs);
-  }
-
-  int ind = 0;
-  for (int fi = 0; fi <numF; ++fi)
-  {
-    const Eigen::Matrix<typename DerivedV::Scalar, 1, 3> &b1 = B1.row(fi);
-    const Eigen::Matrix<typename DerivedV::Scalar, 1, 3> &b2 = B2.row(fi);
-    if(isConstrained[fi])
-    {
-      std::complex<typename DerivedV::Scalar> ck(0);
-
-      for (int j = 0; j < allCombs.rows(); ++j)
-      {
-        std::complex<typename DerivedV::Scalar> tk(1.);
-        //collect products
-        for (int i = 0; i < allCombs.cols(); ++i)
-        {
-          int index = allCombs(j,i);
-
-          int ri = rootsIndex[index];
-          Eigen::Matrix<typename DerivedV::Scalar, 1, 3> w;
-          if (ri>0)
-            w = cfW.block(fi,3*(ri-1),1,3);
-          else
-            w = -cfW.block(fi,3*(-ri-1),1,3);
-          typename DerivedV::Scalar w0 = w.dot(b1);
-          typename DerivedV::Scalar w1 = w.dot(b2);
-          std::complex<typename DerivedV::Scalar> u(w0,w1);
-          tk*= u;
-        }
-        //collect sum
-        ck += tk;
-      }
-      Ck(ind) = ck;
-      ind ++;
-    }
-  }
-
-
-}
-
-template<typename DerivedV, typename DerivedF>
-IGL_INLINE void igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::computek()
-{
-  K.setZero(numE);
-  // For every non-border edge
-  for (unsigned eid=0; eid<numE; ++eid)
-  {
-    if (!isBorderEdge[eid])
-    {
-      int fid0 = E2F(eid,0);
-      int fid1 = E2F(eid,1);
-
-      Eigen::Matrix<typename DerivedV::Scalar, 1, 3> N0 = FN.row(fid0);
-      Eigen::Matrix<typename DerivedV::Scalar, 1, 3> N1 = FN.row(fid1);
-
-      // find common edge on triangle 0 and 1
-      int fid0_vc = -1;
-      int fid1_vc = -1;
-      for (unsigned i=0;i<3;++i)
-      {
-        if (F2E(fid0,i) == eid)
-          fid0_vc = i;
-        if (F2E(fid1,i) == eid)
-          fid1_vc = i;
-      }
-      assert(fid0_vc != -1);
-      assert(fid1_vc != -1);
-
-      Eigen::Matrix<typename DerivedV::Scalar, 1, 3> common_edge = V.row(F(fid0,(fid0_vc+1)%3)) - V.row(F(fid0,fid0_vc));
-      common_edge.normalize();
-
-      // Map the two triangles in a new space where the common edge is the x axis and the N0 the z axis
-      Eigen::Matrix<typename DerivedV::Scalar, 3, 3> P;
-      Eigen::Matrix<typename DerivedV::Scalar, 1, 3> o = V.row(F(fid0,fid0_vc));
-      Eigen::Matrix<typename DerivedV::Scalar, 1, 3> tmp = -N0.cross(common_edge);
-      P << common_edge, tmp, N0;
-//      P.transposeInPlace();
-
-
-      Eigen::Matrix<typename DerivedV::Scalar, 3, 3> V0;
-      V0.row(0) = V.row(F(fid0,0)) -o;
-      V0.row(1) = V.row(F(fid0,1)) -o;
-      V0.row(2) = V.row(F(fid0,2)) -o;
-
-      V0 = (P*V0.transpose()).transpose();
-
-//      assert(V0(0,2) < 1e-10);
-//      assert(V0(1,2) < 1e-10);
-//      assert(V0(2,2) < 1e-10);
-
-      Eigen::Matrix<typename DerivedV::Scalar, 3, 3> V1;
-      V1.row(0) = V.row(F(fid1,0)) -o;
-      V1.row(1) = V.row(F(fid1,1)) -o;
-      V1.row(2) = V.row(F(fid1,2)) -o;
-      V1 = (P*V1.transpose()).transpose();
-
-//      assert(V1(fid1_vc,2) < 10e-10);
-//      assert(V1((fid1_vc+1)%3,2) < 10e-10);
-
-      // compute rotation R such that R * N1 = N0
-      // i.e. map both triangles to the same plane
-      double alpha = -atan2(V1((fid1_vc+2)%3,2),V1((fid1_vc+2)%3,1));
-
-      Eigen::Matrix<typename DerivedV::Scalar, 3, 3> R;
-      R << 1,          0,            0,
-      0, cos(alpha), -sin(alpha) ,
-      0, sin(alpha),  cos(alpha);
-      V1 = (R*V1.transpose()).transpose();
-
-//      assert(V1(0,2) < 1e-10);
-//      assert(V1(1,2) < 1e-10);
-//      assert(V1(2,2) < 1e-10);
-
-      // measure the angle between the reference frames
-      // k_ij is the angle between the triangle on the left and the one on the right
-      Eigen::Matrix<typename DerivedV::Scalar, 1, 3> ref0 = V0.row(1) - V0.row(0);
-      Eigen::Matrix<typename DerivedV::Scalar, 1, 3> ref1 = V1.row(1) - V1.row(0);
-
-      ref0.normalize();
-      ref1.normalize();
-
-      double ktemp = atan2(ref1(1),ref1(0)) - atan2(ref0(1),ref0(0));
-
-      // just to be sure, rotate ref0 using angle ktemp...
-      Eigen::Matrix<typename DerivedV::Scalar, 2, 2> R2;
-      R2 << cos(ktemp), -sin(ktemp), sin(ktemp), cos(ktemp);
-
-      Eigen::Matrix<typename DerivedV::Scalar, 1, 2> tmp1 = R2*(ref0.head(2)).transpose();
-
-//      assert(tmp1(0) - ref1(0) < 1e-10);
-//      assert(tmp1(1) - ref1(1) < 1e-10);
-
-      K[eid] = ktemp;
-    }
-  }
-
-}
+//#include <igl/edge_topology.h>
+//#include <igl/local_basis.h>
+//#include <igl/nchoosek.h>
+//#include <igl/slice.h>
+//#include <igl/polyroots.h>
+//#include <Eigen/Sparse>
+//#include <Eigen/Geometry>
+//#include <iostream>
+//#include <iostream>
+//
+//namespace igl {
+//  template <typename DerivedV, typename DerivedF>
+//  class GeneralPolyVectorFieldFinder
+//  {
+//  private:
+//    const Eigen::PlainObjectBase<DerivedV> &V;
+//    const Eigen::PlainObjectBase<DerivedF> &F; int numF;
+//    const int n;
+//
+//    Eigen::MatrixXi EV; int numE;
+//    Eigen::MatrixXi F2E;
+//    Eigen::MatrixXi E2F;
+//    Eigen::VectorXd K;
+//
+//    Eigen::VectorXi isBorderEdge;
+//    int numInteriorEdges;
+//    Eigen::Matrix<int,Eigen::Dynamic,2> E2F_int;
+//    Eigen::VectorXi indInteriorToFull;
+//    Eigen::VectorXi indFullToInterior;
+//
+////#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
+//    Eigen::PlainObjectBase<DerivedV> B1, B2, FN;
+//
+//    IGL_INLINE void computek();
+//    IGL_INLINE void setFieldFromGeneralCoefficients(const  std::vector<Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1>> &coeffs,
+//                                                    std::vector<Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 2> > &pv);
+//    IGL_INLINE void computeCoefficientLaplacian(int n, Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &D);
+//    IGL_INLINE void getGeneralCoeffConstraints(const Eigen::VectorXi &isConstrained,
+//                                    const Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> &cfW,
+//                                    int k,
+//                                    const Eigen::VectorXi &rootsIndex,
+//                                    Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1> &Ck);
+//    IGL_INLINE void precomputeInteriorEdges();
+//
+//
+//    IGL_INLINE void minQuadWithKnownMini(const Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &Q,
+//                                         const Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &f,
+//                                         const Eigen::VectorXi isConstrained,
+//                                         const Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> &xknown,
+//                                         Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> &x);
+//
+//  public:
+//    IGL_INLINE GeneralPolyVectorFieldFinder(const Eigen::PlainObjectBase<DerivedV> &_V,
+//                                     const Eigen::PlainObjectBase<DerivedF> &_F,
+//                                     const int &_n);
+//    IGL_INLINE bool solve(const Eigen::VectorXi &isConstrained,
+//               const Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> &cfW,
+//               const Eigen::VectorXi &rootsIndex,
+//               Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> &output);
+//
+//  };
+//}
+//
+//template<typename DerivedV, typename DerivedF>
+//IGL_INLINE igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::
+//          GeneralPolyVectorFieldFinder(const Eigen::PlainObjectBase<DerivedV> &_V,
+//                                const Eigen::PlainObjectBase<DerivedF> &_F,
+//                                const int &_n):
+//V(_V),
+//F(_F),
+//numF(_F.rows()),
+//n(_n)
+//{
+//
+//  igl::edge_topology(V,F,EV,F2E,E2F);
+//  numE = EV.rows();
+//
+//
+//  precomputeInteriorEdges();
+//
+//  igl::local_basis(V,F,B1,B2,FN);
+//
+//  computek();
+//
+//};
+//
+//
+//template<typename DerivedV, typename DerivedF>
+//IGL_INLINE void igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::
+//precomputeInteriorEdges()
+//{
+//  // Flag border edges
+//  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)))
+//      isBorderEdge[i] = 1;
+//      else
+//      {
+//        indFullToInterior[i] = numInteriorEdges;
+//        numInteriorEdges++;
+//      }
+//  }
+//
+//  E2F_int.resize(numInteriorEdges, 2);
+//  indInteriorToFull.setZero(numInteriorEdges,1);
+//  int ii = 0;
+//  for (int k=0; k<numE; ++k)
+//  {
+//    if (isBorderEdge[k])
+//      continue;
+//    E2F_int.row(ii) = E2F.row(k);
+//    indInteriorToFull[ii] = k;
+//    ii++;
+//  }
+//
+//}
+//
+//
+//template<typename DerivedV, typename DerivedF>
+//IGL_INLINE void igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::
+//minQuadWithKnownMini(const Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &Q,
+//                          const Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &f,
+//                     const Eigen::VectorXi isConstrained,
+//                          const Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> &xknown,
+//                          Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> &x)
+//{
+//  int N = Q.rows();
+//
+//  int nc = xknown.rows();
+//  Eigen::VectorXi known; known.setZero(nc,1);
+//  Eigen::VectorXi unknown; unknown.setZero(N-nc,1);
+//
+//  int indk = 0, indu = 0;
+//  for (int i = 0; i<N; ++i)
+//    if (isConstrained[i])
+//    {
+//      known[indk] = i;
+//      indk++;
+//    }
+//    else
+//    {
+//      unknown[indu] = i;
+//      indu++;
+//    }
+//
+//  Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar>> Quu, Quk;
+//
+//  igl::slice(Q,unknown, unknown, Quu);
+//  igl::slice(Q,unknown, known, Quk);
+//
+//
+//  std::vector<typename Eigen::Triplet<std::complex<typename DerivedV::Scalar> > > tripletList;
+//
+//  Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > fu(N-nc,1);
+//
+//  igl::slice(f,unknown, Eigen::VectorXi::Zero(1,1), fu);
+//
+//  Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > rhs = (Quk*xknown).sparseView()+.5*fu;
+//
+//  Eigen::SparseLU< Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar>>> solver;
+//  solver.compute(-Quu);
+//  if(solver.info()!=Eigen::Success)
+//  {
+//    std::cerr<<"Decomposition failed!"<<std::endl;
+//    return;
+//  }
+//  Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar>>  b  = solver.solve(rhs);
+//  if(solver.info()!=Eigen::Success)
+//  {
+//    std::cerr<<"Solving failed!"<<std::endl;
+//    return;
+//  }
+//
+//  indk = 0, indu = 0;
+//  x.setZero(N,1);
+//  for (int i = 0; i<N; ++i)
+//    if (isConstrained[i])
+//      x[i] = xknown[indk++];
+//    else
+//      x[i] = b.coeff(indu++,0);
+//
+//}
+//
+//
+//
+//template<typename DerivedV, typename DerivedF>
+//IGL_INLINE bool igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::
+//                     solve(const Eigen::VectorXi &isConstrained,
+//                           const Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> &cfW,
+//                           const Eigen::VectorXi &rootsIndex,
+//                           Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> &output)
+//{
+//	std::cerr << "This code is broken!" << std::endl;
+//	exit(1);
+//  // polynomial is of the form:
+//  // z^(2n) +
+//  // -c[0]z^(2n-1) +
+//  // c[1]z^(2n-2) +
+//  // -c[2]z^(2n-3) +
+//  // ... +
+//  // (-1)^n c[n-1]
+//  //std::vector<Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1>> coeffs(n,Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1>::Zero(numF, 1));
+//  //for (int i =0; i<n; ++i)
+//  //{
+//  //  int degree = i+1;
+//  //  Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1> Ck;
+//  //  getGeneralCoeffConstraints(isConstrained,
+//  //                             cfW,
+//  //                             i,
+//  //                             rootsIndex,
+//  //                             Ck);
+//
+//  //  Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > DD;
+//  //  computeCoefficientLaplacian(degree, DD);
+//  //  Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > f; f.resize(numF,1);
+//
+//  //  if (isConstrained.sum() == numF)
+//  //    coeffs[i] = Ck;
+//  //  else
+//  //    minQuadWithKnownMini(DD, f, isConstrained, Ck, coeffs[i]);
+//  //}
+//  //std::vector<Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 2> > pv;
+//  //setFieldFromGeneralCoefficients(coeffs, pv);
+//  //output.setZero(numF,3*n);
+//  //for (int fi=0; fi<numF; ++fi)
+//  //{
+//  //  const Eigen::Matrix<typename DerivedV::Scalar, 1, 3> &b1 = B1.row(fi);
+//  //  const Eigen::Matrix<typename DerivedV::Scalar, 1, 3> &b2 = B2.row(fi);
+//  //  for (int i=0; i<n; ++i)
+//  //    output.block(fi,3*i, 1, 3) = pv[i](fi,0)*b1 + pv[i](fi,1)*b2;
+//  //}
+//  return true;
+//}
+//
+//template<typename DerivedV, typename DerivedF>
+//IGL_INLINE void igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::setFieldFromGeneralCoefficients(const  std::vector<Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1>> &coeffs,
+//                                                            std::vector<Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 2>> &pv)
+//{
+//  pv.assign(n, Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 2>::Zero(numF, 2));
+//  for (int i = 0; i <numF; ++i)
+//  {
+//
+//    //    poly coefficients: 1, 0, -Acoeff, 0, Bcoeff
+//    //    matlab code from roots (given there are no trailing zeros in the polynomial coefficients)
+//    Eigen::Matrix<std::complex<DerivedV::Scalar>, Eigen::Dynamic,1> polyCoeff;
+//    polyCoeff.setZero(n+1,1);
+//    polyCoeff[0] = 1.;
+//    int sign = 1;
+//    for (int k =0; k<n; ++k)
+//    {
+//      sign = -sign;
+//      int degree = k+1;
+//      polyCoeff[degree] = (1.*sign)*coeffs[k](i);
+//    }
+//
+//    Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1> roots;
+//    igl::polyRoots<std::complex<typename DerivedV::Scalar>, typename DerivedV::Scalar >(polyCoeff,roots);
+//    for (int k=0; k<n; ++k)
+//    {
+//      pv[k](i,0) = real(roots[k]);
+//      pv[k](i,1) = imag(roots[k]);
+//    }
+//  }
+//
+//}
+//
+//
+//template<typename DerivedV, typename DerivedF>
+//IGL_INLINE void igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::computeCoefficientLaplacian(int n, Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &D)
+//{
+//  std::vector<Eigen::Triplet<std::complex<typename DerivedV::Scalar> >> tripletList;
+//
+//  // For every non-border edge
+//  for (unsigned eid=0; eid<numE; ++eid)
+//  {
+//    if (!isBorderEdge[eid])
+//    {
+//      int fid0 = E2F(eid,0);
+//      int fid1 = E2F(eid,1);
+//
+//      tripletList.push_back(Eigen::Triplet<std::complex<typename DerivedV::Scalar> >(fid0,
+//                                           fid0,
+//                                           std::complex<typename DerivedV::Scalar>(1.)));
+//      tripletList.push_back(Eigen::Triplet<std::complex<typename DerivedV::Scalar> >(fid1,
+//                                           fid1,
+//                                           std::complex<typename DerivedV::Scalar>(1.)));
+//      tripletList.push_back(Eigen::Triplet<std::complex<typename DerivedV::Scalar> >(fid0,
+//                                           fid1,
+//                                                                                     -1.*std::polar(1.,-1.*n*K[eid])));
+//      tripletList.push_back(Eigen::Triplet<std::complex<typename DerivedV::Scalar> >(fid1,
+//                                           fid0,
+//                                                                                     -1.*std::polar(1.,1.*n*K[eid])));
+//
+//    }
+//  }
+//  D.resize(numF,numF);
+//  D.setFromTriplets(tripletList.begin(), tripletList.end());
+//
+//
+//}
+//
+////this gives the coefficients without the (-1)^k that multiplies them
+//template<typename DerivedV, typename DerivedF>
+//IGL_INLINE void igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::getGeneralCoeffConstraints(const Eigen::VectorXi &isConstrained,
+//                                                       const Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> &cfW,
+//                                                       int k,
+//                                                       const Eigen::VectorXi &rootsIndex,
+//                                                       Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1> &Ck)
+//{
+//  int numConstrained = isConstrained.sum();
+//  Ck.resize(numConstrained,1);
+//  // int n = rootsIndex.cols();
+//
+//  Eigen::MatrixXi allCombs;
+//  {
+//    Eigen::VectorXi V = Eigen::VectorXi::LinSpaced(n,0,n-1);
+//    igl::nchoosek(V,k+1,allCombs);
+//  }
+//
+//  int ind = 0;
+//  for (int fi = 0; fi <numF; ++fi)
+//  {
+//    const Eigen::Matrix<typename DerivedV::Scalar, 1, 3> &b1 = B1.row(fi);
+//    const Eigen::Matrix<typename DerivedV::Scalar, 1, 3> &b2 = B2.row(fi);
+//    if(isConstrained[fi])
+//    {
+//      std::complex<typename DerivedV::Scalar> ck(0);
+//
+//      for (int j = 0; j < allCombs.rows(); ++j)
+//      {
+//        std::complex<typename DerivedV::Scalar> tk(1.);
+//        //collect products
+//        for (int i = 0; i < allCombs.cols(); ++i)
+//        {
+//          int index = allCombs(j,i);
+//
+//          int ri = rootsIndex[index];
+//          Eigen::Matrix<typename DerivedV::Scalar, 1, 3> w;
+//          if (ri>0)
+//            w = cfW.block(fi,3*(ri-1),1,3);
+//          else
+//            w = -cfW.block(fi,3*(-ri-1),1,3);
+//          typename DerivedV::Scalar w0 = w.dot(b1);
+//          typename DerivedV::Scalar w1 = w.dot(b2);
+//          std::complex<typename DerivedV::Scalar> u(w0,w1);
+//          tk*= u;
+//        }
+//        //collect sum
+//        ck += tk;
+//      }
+//      Ck(ind) = ck;
+//      ind ++;
+//    }
+//  }
+//
+//
+//}
+//
+//template<typename DerivedV, typename DerivedF>
+//IGL_INLINE void igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::computek()
+//{
+//  K.setZero(numE);
+//  // For every non-border edge
+//  for (unsigned eid=0; eid<numE; ++eid)
+//  {
+//    if (!isBorderEdge[eid])
+//    {
+//      int fid0 = E2F(eid,0);
+//      int fid1 = E2F(eid,1);
+//
+//      Eigen::Matrix<typename DerivedV::Scalar, 1, 3> N0 = FN.row(fid0);
+//      Eigen::Matrix<typename DerivedV::Scalar, 1, 3> N1 = FN.row(fid1);
+//
+//      // find common edge on triangle 0 and 1
+//      int fid0_vc = -1;
+//      int fid1_vc = -1;
+//      for (unsigned i=0;i<3;++i)
+//      {
+//        if (F2E(fid0,i) == eid)
+//          fid0_vc = i;
+//        if (F2E(fid1,i) == eid)
+//          fid1_vc = i;
+//      }
+//      assert(fid0_vc != -1);
+//      assert(fid1_vc != -1);
+//
+//      Eigen::Matrix<typename DerivedV::Scalar, 1, 3> common_edge = V.row(F(fid0,(fid0_vc+1)%3)) - V.row(F(fid0,fid0_vc));
+//      common_edge.normalize();
+//
+//      // Map the two triangles in a new space where the common edge is the x axis and the N0 the z axis
+//      Eigen::Matrix<typename DerivedV::Scalar, 3, 3> P;
+//      Eigen::Matrix<typename DerivedV::Scalar, 1, 3> o = V.row(F(fid0,fid0_vc));
+//      Eigen::Matrix<typename DerivedV::Scalar, 1, 3> tmp = -N0.cross(common_edge);
+//      P << common_edge, tmp, N0;
+////      P.transposeInPlace();
+//
+//
+//      Eigen::Matrix<typename DerivedV::Scalar, 3, 3> V0;
+//      V0.row(0) = V.row(F(fid0,0)) -o;
+//      V0.row(1) = V.row(F(fid0,1)) -o;
+//      V0.row(2) = V.row(F(fid0,2)) -o;
+//
+//      V0 = (P*V0.transpose()).transpose();
+//
+////      assert(V0(0,2) < 1e-10);
+////      assert(V0(1,2) < 1e-10);
+////      assert(V0(2,2) < 1e-10);
+//
+//      Eigen::Matrix<typename DerivedV::Scalar, 3, 3> V1;
+//      V1.row(0) = V.row(F(fid1,0)) -o;
+//      V1.row(1) = V.row(F(fid1,1)) -o;
+//      V1.row(2) = V.row(F(fid1,2)) -o;
+//      V1 = (P*V1.transpose()).transpose();
+//
+////      assert(V1(fid1_vc,2) < 10e-10);
+////      assert(V1((fid1_vc+1)%3,2) < 10e-10);
+//
+//      // compute rotation R such that R * N1 = N0
+//      // i.e. map both triangles to the same plane
+//      double alpha = -atan2(V1((fid1_vc+2)%3,2),V1((fid1_vc+2)%3,1));
+//
+//      Eigen::Matrix<typename DerivedV::Scalar, 3, 3> R;
+//      R << 1,          0,            0,
+//      0, cos(alpha), -sin(alpha) ,
+//      0, sin(alpha),  cos(alpha);
+//      V1 = (R*V1.transpose()).transpose();
+//
+////      assert(V1(0,2) < 1e-10);
+////      assert(V1(1,2) < 1e-10);
+////      assert(V1(2,2) < 1e-10);
+//
+//      // measure the angle between the reference frames
+//      // k_ij is the angle between the triangle on the left and the one on the right
+//      Eigen::Matrix<typename DerivedV::Scalar, 1, 3> ref0 = V0.row(1) - V0.row(0);
+//      Eigen::Matrix<typename DerivedV::Scalar, 1, 3> ref1 = V1.row(1) - V1.row(0);
+//
+//      ref0.normalize();
+//      ref1.normalize();
+//
+//      double ktemp = atan2(ref1(1),ref1(0)) - atan2(ref0(1),ref0(0));
+//
+//      // just to be sure, rotate ref0 using angle ktemp...
+//      Eigen::Matrix<typename DerivedV::Scalar, 2, 2> R2;
+//      R2 << cos(ktemp), -sin(ktemp), sin(ktemp), cos(ktemp);
+//
+//      Eigen::Matrix<typename DerivedV::Scalar, 1, 2> tmp1 = R2*(ref0.head(2)).transpose();
+//
+////      assert(tmp1(0) - ref1(0) < 1e-10);
+////      assert(tmp1(1) - ref1(1) < 1e-10);
+//
+//      K[eid] = ktemp;
+//    }
+//  }
+//
+//}
 
 
 IGL_INLINE void igl::n_polyvector_general(const Eigen::MatrixXd &V,
@@ -478,17 +475,19 @@ IGL_INLINE void igl::n_polyvector_general(const Eigen::MatrixXd &V,
                              const Eigen::VectorXi &I,
                              Eigen::MatrixXd &output)
 {
-  Eigen::VectorXi isConstrained = Eigen::VectorXi::Constant(F.rows(),0);
-  Eigen::MatrixXd cfW = Eigen::MatrixXd::Constant(F.rows(),bc.cols(),0);
-
-  for(unsigned i=0; i<b.size();++i)
-  {
-    isConstrained(b(i)) = 1;
-    cfW.row(b(i)) << bc.row(i);
-  }
-  int n = I.rows();
-  igl::GeneralPolyVectorFieldFinder<Eigen::MatrixXd, Eigen::MatrixXi> pvff(V,F,n);
-  pvff.solve(isConstrained, cfW, I, output);
+	// This functions is broken, please contact Olga Diamanti
+	assert(0);
+  //Eigen::VectorXi isConstrained = Eigen::VectorXi::Constant(F.rows(),0);
+  //Eigen::MatrixXd cfW = Eigen::MatrixXd::Constant(F.rows(),bc.cols(),0);
+
+  //for(unsigned i=0; i<b.size();++i)
+  //{
+  //  isConstrained(b(i)) = 1;
+  //  cfW.row(b(i)) << bc.row(i);
+  //}
+  //int n = I.rows();
+  //igl::GeneralPolyVectorFieldFinder<Eigen::MatrixXd, Eigen::MatrixXi> pvff(V,F,n);
+  //pvff.solve(isConstrained, cfW, I, output);
 }
 
 

include/igl/polyvector_field_matchings.cpp

@@ -168,7 +168,7 @@ IGL_INLINE typename DerivedC::Scalar igl::polyvector_field_matchings(
   Eigen::MatrixXi E, E2F, F2E;
   igl::edge_topology(V,F,E,F2E,E2F);
 
-#warning "Poor templating of igl::polyvector_field_matchings forces FN to be DerivedV (this could cause issues if DerivedV has fixed number of rows)"
+//#warning "Poor templating of igl::polyvector_field_matchings forces FN to be DerivedV (this could cause issues if DerivedV has fixed number of rows)"
   DerivedV FN;
   igl::per_face_normals(V,F,FN);