Merge branch 'master' into bugfix-get_smallest_pos_quad_zero

Former-commit-id: 896cf7a23e3cbebcbdd4b6b3243ee145a7283bd2

include/igl/AABB.cpp

@@ -148,22 +148,14 @@ IGL_INLINE void igl::AABB<DerivedV,DIM>::init(
         {
           SIdI(i) = SI(I(i),max_d);
         }
-        // Since later I use <= I think I don't need to worry about odd/even
         // Pass by copy to avoid changing input
-        const auto median = [](VectorXi A)->Scalar
+        const auto median = [](VectorXi A)->int
         {
-          size_t n = A.size()/2;
+          size_t n = (A.size()-1)/2;
           nth_element(A.data(),A.data()+n,A.data()+A.size());
-          if(A.rows() % 2 == 1)
-          {
-            return A(n);
-          }else
-          {
-            nth_element(A.data(),A.data()+n-1,A.data()+A.size());
-            return 0.5*(A(n)+A(n-1));
-          }
+          return A(n);
         };
-        const Scalar med = median(SIdI);
+        const int med = median(SIdI);
         VectorXi LI((I.rows()+1)/2),RI(I.rows()/2);
         assert(LI.rows()+RI.rows() == I.rows());
         // Distribute left and right

include/igl/HalfEdgeIterator.cpp

@@ -8,11 +8,11 @@
 
 #include "HalfEdgeIterator.h"
 
-template <typename DerivedF>
-IGL_INLINE igl::HalfEdgeIterator<DerivedF>::HalfEdgeIterator(
+template <typename DerivedF, typename DerivedFF, typename DerivedFFi>
+IGL_INLINE igl::HalfEdgeIterator<DerivedF,DerivedFF,DerivedFFi>::HalfEdgeIterator(
     const Eigen::PlainObjectBase<DerivedF>& _F,
-    const Eigen::PlainObjectBase<DerivedF>& _FF,
-    const Eigen::PlainObjectBase<DerivedF>& _FFi,
+    const Eigen::PlainObjectBase<DerivedFF>& _FF,
+    const Eigen::PlainObjectBase<DerivedFFi>& _FFi,
     int _fi,
     int _ei,
     bool _reverse
@@ -20,8 +20,8 @@ IGL_INLINE igl::HalfEdgeIterator<DerivedF>::HalfEdgeIterator(
 : fi(_fi), ei(_ei), reverse(_reverse), F(_F), FF(_FF), FFi(_FFi)
 {}
 
-template <typename DerivedF>
-IGL_INLINE void igl::HalfEdgeIterator<DerivedF>::flipF()
+template <typename DerivedF, typename DerivedFF, typename DerivedFFi>
+IGL_INLINE void igl::HalfEdgeIterator<DerivedF,DerivedFF,DerivedFFi>::flipF()
 {
   if (isBorder())
     return;
@@ -36,8 +36,8 @@ IGL_INLINE void igl::HalfEdgeIterator<DerivedF>::flipF()
 
 
 // Change Edge
-template <typename DerivedF>
-IGL_INLINE void igl::HalfEdgeIterator<DerivedF>::flipE()
+template <typename DerivedF, typename DerivedFF, typename DerivedFFi>
+IGL_INLINE void igl::HalfEdgeIterator<DerivedF,DerivedFF,DerivedFFi>::flipE()
 {
   if (!reverse)
     ei = (ei+2)%3; // ei-1
@@ -48,14 +48,14 @@ IGL_INLINE void igl::HalfEdgeIterator<DerivedF>::flipE()
 }
 
 // Change Vertex
-template <typename DerivedF>
-IGL_INLINE void igl::HalfEdgeIterator<DerivedF>::flipV()
+template <typename DerivedF, typename DerivedFF, typename DerivedFFi>
+IGL_INLINE void igl::HalfEdgeIterator<DerivedF,DerivedFF,DerivedFFi>::flipV()
 {
   reverse = !reverse;
 }
 
-template <typename DerivedF>
-IGL_INLINE bool igl::HalfEdgeIterator<DerivedF>::isBorder()
+template <typename DerivedF, typename DerivedFF, typename DerivedFFi>
+IGL_INLINE bool igl::HalfEdgeIterator<DerivedF,DerivedFF,DerivedFFi>::isBorder()
 {
   return (FF)(fi,ei) == -1;
 }
@@ -71,8 +71,8 @@ IGL_INLINE bool igl::HalfEdgeIterator<DerivedF>::isBorder()
  * In this example, if a and d are of-border and the pos is iterating counterclockwise, this method iterate through the faces incident on vertex v,
  * producing the sequence a, b, c, d, a, b, c, ...
  */
-template <typename DerivedF>
-IGL_INLINE bool igl::HalfEdgeIterator<DerivedF>::NextFE()
+template <typename DerivedF, typename DerivedFF, typename DerivedFFi>
+IGL_INLINE bool igl::HalfEdgeIterator<DerivedF,DerivedFF,DerivedFFi>::NextFE()
 {
   if ( isBorder() ) // we are on a border
   {
@@ -93,8 +93,8 @@ IGL_INLINE bool igl::HalfEdgeIterator<DerivedF>::NextFE()
 }
 
 // Get vertex index
-template <typename DerivedF>
-IGL_INLINE int igl::HalfEdgeIterator<DerivedF>::Vi()
+template <typename DerivedF, typename DerivedFF, typename DerivedFFi>
+IGL_INLINE int igl::HalfEdgeIterator<DerivedF,DerivedFF,DerivedFFi>::Vi()
 {
   assert(fi >= 0);
   assert(fi < F.rows());
@@ -108,22 +108,22 @@ IGL_INLINE int igl::HalfEdgeIterator<DerivedF>::Vi()
 }
 
 // Get face index
-template <typename DerivedF>
-IGL_INLINE int igl::HalfEdgeIterator<DerivedF>::Fi()
+template <typename DerivedF, typename DerivedFF, typename DerivedFFi>
+IGL_INLINE int igl::HalfEdgeIterator<DerivedF,DerivedFF,DerivedFFi>::Fi()
 {
   return fi;
 }
 
 // Get edge index
-template <typename DerivedF>
-IGL_INLINE int igl::HalfEdgeIterator<DerivedF>::Ei()
+template <typename DerivedF, typename DerivedFF, typename DerivedFFi>
+IGL_INLINE int igl::HalfEdgeIterator<DerivedF,DerivedFF,DerivedFFi>::Ei()
 {
   return ei;
 }
 
 
-template <typename DerivedF>
-IGL_INLINE bool igl::HalfEdgeIterator<DerivedF>::operator==(HalfEdgeIterator& p2)
+template <typename DerivedF, typename DerivedFF, typename DerivedFFi>
+IGL_INLINE bool igl::HalfEdgeIterator<DerivedF,DerivedFF,DerivedFFi>::operator==(HalfEdgeIterator& p2)
 {
   return
       (
@@ -138,17 +138,21 @@ IGL_INLINE bool igl::HalfEdgeIterator<DerivedF>::operator==(HalfEdgeIterator& p2
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
-template igl::HalfEdgeIterator<Eigen::Matrix<int, -1, 3, 0, -1, 3> >::HalfEdgeIterator(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, int, int, bool);
-template bool igl::HalfEdgeIterator<Eigen::Matrix<int, -1, -1, 0, -1, -1> >::NextFE();
-template int igl::HalfEdgeIterator<Eigen::Matrix<int, -1, 3, 0, -1, 3> >::Ei();
-template int igl::HalfEdgeIterator<Eigen::Matrix<int, -1, 3, 0, -1, 3> >::Fi();
-template bool igl::HalfEdgeIterator<Eigen::Matrix<int, -1, 3, 0, -1, 3> >::NextFE();
-template int igl::HalfEdgeIterator<Eigen::Matrix<int, -1, -1, 0, -1, -1> >::Vi();
-template igl::HalfEdgeIterator<Eigen::Matrix<int, -1, -1, 0, -1, -1> >::HalfEdgeIterator(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int, int, bool);
-template int igl::HalfEdgeIterator<Eigen::Matrix<int, -1, -1, 0, -1, -1> >::Ei();
-template int igl::HalfEdgeIterator<Eigen::Matrix<int, -1, -1, 0, -1, -1> >::Fi();
-template void igl::HalfEdgeIterator<Eigen::Matrix<int, -1, -1, 0, -1, -1> >::flipE();
-template void igl::HalfEdgeIterator<Eigen::Matrix<int, -1, -1, 0, -1, -1> >::flipF();
-template void igl::HalfEdgeIterator<Eigen::Matrix<int, -1, -1, 0, -1, -1> >::flipV();
-template bool igl::HalfEdgeIterator<Eigen::Matrix<int, -1, -1, 0, -1, -1> >::operator==(igl::HalfEdgeIterator<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
-#endif
+template      igl::HalfEdgeIterator<Eigen::Matrix<int, -1, 3, 0, -1, 3>  ,Eigen::Matrix<int, -1, 3, 0, -1, 3>  ,Eigen::Matrix<int, -1, 3, 0, -1, 3>   >::HalfEdgeIterator(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, int, int, bool);
+template igl::HalfEdgeIterator<Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >::HalfEdgeIterator(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int, int, bool);
+template bool igl::HalfEdgeIterator<Eigen::Matrix<int, -1, -1, 0, -1, -1>,Eigen::Matrix<int, -1, -1, 0, -1, -1>,Eigen::Matrix<int, -1, -1, 0, -1, -1> >::NextFE();
+template int  igl::HalfEdgeIterator<Eigen::Matrix<int, -1, -1, 0, -1, -1>,Eigen::Matrix<int, -1, -1, 0, -1, -1>,Eigen::Matrix<int, -1, -1, 0, -1, -1> >::Ei();
+template int  igl::HalfEdgeIterator<Eigen::Matrix<int, -1, 3, 0, -1, 3>  ,Eigen::Matrix<int, -1, -1, 0, -1, -1>,Eigen::Matrix<int, -1, -1, 0, -1, -1> >::Ei();
+template int  igl::HalfEdgeIterator<Eigen::Matrix<int, -1, 3, 0, -1, 3>  ,Eigen::Matrix<int, -1, 3, 0, -1, 3>  ,Eigen::Matrix<int, -1, 3, 0, -1, 3>   >::Ei();
+template int  igl::HalfEdgeIterator<Eigen::Matrix<int, -1, 3, 0, -1, 3>  ,Eigen::Matrix<int, -1, 3, 0, -1, 3>  ,Eigen::Matrix<int, -1, 3, 0, -1, 3>   >::Fi();
+template bool igl::HalfEdgeIterator<Eigen::Matrix<int, -1, 3, 0, -1, 3>  ,Eigen::Matrix<int, -1, 3, 0, -1, 3>  ,Eigen::Matrix<int, -1, 3, 0, -1, 3>   >::NextFE();
+template int  igl::HalfEdgeIterator<Eigen::Matrix<int, -1, -1, 0, -1, -1>,Eigen::Matrix<int, -1, -1, 0, -1, -1>,Eigen::Matrix<int, -1, -1, 0, -1, -1> >::Vi();
+template      igl::HalfEdgeIterator<Eigen::Matrix<int, -1, -1, 0, -1, -1>,Eigen::Matrix<int, -1, -1, 0, -1, -1>,Eigen::Matrix<int, -1, -1, 0, -1, -1> >::HalfEdgeIterator(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int, int, bool);
+template int  igl::HalfEdgeIterator<Eigen::Matrix<int, -1, -1, 0, -1, -1>,Eigen::Matrix<int, -1, -1, 0, -1, -1>,Eigen::Matrix<int, -1, -1, 0, -1, -1> >::Fi();
+template void igl::HalfEdgeIterator<Eigen::Matrix<int, -1, -1, 0, -1, -1>,Eigen::Matrix<int, -1, -1, 0, -1, -1>,Eigen::Matrix<int, -1, -1, 0, -1, -1> >::flipE();
+template void igl::HalfEdgeIterator<Eigen::Matrix<int, -1, -1, 0, -1, -1>,Eigen::Matrix<int, -1, -1, 0, -1, -1>,Eigen::Matrix<int, -1, -1, 0, -1, -1> >::flipF();
+template void igl::HalfEdgeIterator<Eigen::Matrix<int, -1, -1, 0, -1, -1>,Eigen::Matrix<int, -1, -1, 0, -1, -1>,Eigen::Matrix<int, -1, -1, 0, -1, -1> >::flipV();
+template bool igl::HalfEdgeIterator<Eigen::Matrix<int, -1, -1, 0, -1, -1>,Eigen::Matrix<int, -1, -1, 0, -1, -1>,Eigen::Matrix<int, -1, -1, 0, -1, -1> >::operator==(igl::HalfEdgeIterator<Eigen::Matrix<int, -1, -1, 0, -1, -1>,Eigen::Matrix<int, -1, -1, 0, -1, -1>,Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template int igl::HalfEdgeIterator<Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >::Fi();
+template bool igl::HalfEdgeIterator<Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >::NextFE();
+#endif

include/igl/HalfEdgeIterator.h

@@ -13,6 +13,8 @@
 #include <vector>
 #include <igl/igl_inline.h>
 
+// This file violates many of the libigl style guidelines.
+
 namespace igl
 {
   // HalfEdgeIterator - Fake halfedge for fast and easy navigation
@@ -29,8 +31,6 @@ namespace igl
   // Each tuple contains information on (face, edge, vertex)
   //    and encoded by (face, edge \in {0,1,2}, bool reverse)
   //
-  // Templates:
-  //    DerivedF Matrix Type for F. Has to be explicitly declared.
   // Inputs:
   //    F #F by 3 list of "faces"
   //    FF #F by 3 list of triangle-triangle adjacency.
@@ -40,15 +40,18 @@ namespace igl
   //    FlipF/E/V changes solely one actual face/edge/vertex resp.
   //    NextFE iterates through one-ring of a vertex robustly.
   //
-  template <typename DerivedF>
+  template <
+    typename DerivedF,
+    typename DerivedFF,
+    typename DerivedFFi>
   class HalfEdgeIterator
   {
   public:
     // Init the HalfEdgeIterator by specifying Face,Edge Index and Orientation
     IGL_INLINE HalfEdgeIterator(
         const Eigen::PlainObjectBase<DerivedF>& _F,
-        const Eigen::PlainObjectBase<DerivedF>& _FF,
-        const Eigen::PlainObjectBase<DerivedF>& _FFi,
+        const Eigen::PlainObjectBase<DerivedFF>& _FF,
+        const Eigen::PlainObjectBase<DerivedFFi>& _FFi,
         int _fi,
         int _ei,
         bool _reverse = false
@@ -73,7 +76,9 @@ namespace igl
      *   / d  \ | / a  \
      *  /______\|/______\
      *          v
-     * In this example, if a and d are of-border and the pos is iterating counterclockwise, this method iterate through the faces incident on vertex v,
+     * In this example, if a and d are of-border and the pos is iterating
+     counterclockwise, this method iterate through the faces incident on vertex
+     v,
      * producing the sequence a, b, c, d, a, b, c, ...
      */
     IGL_INLINE bool NextFE();
@@ -94,9 +99,10 @@ namespace igl
     int ei;
     bool reverse;
 
-    const Eigen::PlainObjectBase<DerivedF>& F;
-    const Eigen::PlainObjectBase<DerivedF>& FF;
-    const Eigen::PlainObjectBase<DerivedF>& FFi;
+    // All the same type? This is likely to break.
+    const DerivedF & F;
+    const DerivedFF & FF;
+    const DerivedFFi & FFi;
   };
 
 }

include/igl/LinSpaced.h

@@ -29,7 +29,8 @@
 namespace igl
 {
   template <typename Derived>
-  inline typename Eigen::DenseBase< Derived >::RandomAccessLinSpacedReturnType LinSpaced(
+  //inline typename Eigen::DenseBase< Derived >::RandomAccessLinSpacedReturnType 
+  inline Derived LinSpaced(
     typename Derived::Index size,
     const typename Derived::Scalar & low,
     const typename Derived::Scalar & high);
@@ -38,7 +39,8 @@ namespace igl
 // Implementation
 
 template <typename Derived>
-inline typename Eigen::DenseBase< Derived >::RandomAccessLinSpacedReturnType 
+//inline typename Eigen::DenseBase< Derived >::RandomAccessLinSpacedReturnType 
+inline Derived
 igl::LinSpaced(
   typename Derived::Index size,
   const typename Derived::Scalar & low,
@@ -48,8 +50,10 @@ igl::LinSpaced(
   {
     // Force empty vector with correct "RandomAccessLinSpacedReturnType" type.
     return Derived::LinSpaced(0,0,1);
-  }else
+  }else if(high < low)
   {
+    return low-Derived::LinSpaced(size,low-low,low-high).array();
+  }else{
     return Derived::LinSpaced(size,low,high);
   }
 }

include/igl/active_set.cpp

@@ -170,10 +170,10 @@ IGL_INLINE igl::SolverStatus igl::active_set(
       old_Z = Z;
     }
 
-    const int as_lx_count = count(as_lx.data(),as_lx.data()+n,TRUE);
-    const int as_ux_count = count(as_ux.data(),as_ux.data()+n,TRUE);
+    const int as_lx_count = std::count(as_lx.data(),as_lx.data()+n,TRUE);
+    const int as_ux_count = std::count(as_ux.data(),as_ux.data()+n,TRUE);
     const int as_ieq_count =
-      count(as_ieq.data(),as_ieq.data()+as_ieq.size(),TRUE);
+      std::count(as_ieq.data(),as_ieq.data()+as_ieq.size(),TRUE);
 #ifndef NDEBUG
     {
       int count = 0;

include/igl/all_edges.cpp

@@ -6,49 +6,21 @@
 // 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 "all_edges.h"
+#include "oriented_facets.h"
 
 template <typename DerivedF, typename DerivedE>
 IGL_INLINE void igl::all_edges(
-  const Eigen::PlainObjectBase<DerivedF> & F,
+  const Eigen::MatrixBase<DerivedF> & F,
   Eigen::PlainObjectBase<DerivedE> & E)
 {
-  E.resize(F.rows()*F.cols(),F.cols()-1);
-  typedef typename DerivedE::Scalar EScalar;
-  switch(F.cols())
-  {
-    case 4:
-      E.block(0*F.rows(),0,F.rows(),1) = F.col(1).template cast<EScalar>();
-      E.block(0*F.rows(),1,F.rows(),1) = F.col(3).template cast<EScalar>();
-      E.block(0*F.rows(),2,F.rows(),1) = F.col(2).template cast<EScalar>();
-
-      E.block(1*F.rows(),0,F.rows(),1) = F.col(0).template cast<EScalar>();
-      E.block(1*F.rows(),1,F.rows(),1) = F.col(2).template cast<EScalar>();
-      E.block(1*F.rows(),2,F.rows(),1) = F.col(3).template cast<EScalar>();
-
-      E.block(2*F.rows(),0,F.rows(),1) = F.col(0).template cast<EScalar>();
-      E.block(2*F.rows(),1,F.rows(),1) = F.col(3).template cast<EScalar>();
-      E.block(2*F.rows(),2,F.rows(),1) = F.col(1).template cast<EScalar>();
-
-      E.block(3*F.rows(),0,F.rows(),1) = F.col(0).template cast<EScalar>();
-      E.block(3*F.rows(),1,F.rows(),1) = F.col(1).template cast<EScalar>();
-      E.block(3*F.rows(),2,F.rows(),1) = F.col(2).template cast<EScalar>();
-      return;
-    case 3:
-      E.block(0*F.rows(),0,F.rows(),1) = F.col(1).template cast<EScalar>();
-      E.block(0*F.rows(),1,F.rows(),1) = F.col(2).template cast<EScalar>();
-      E.block(1*F.rows(),0,F.rows(),1) = F.col(2).template cast<EScalar>();
-      E.block(1*F.rows(),1,F.rows(),1) = F.col(0).template cast<EScalar>();
-      E.block(2*F.rows(),0,F.rows(),1) = F.col(0).template cast<EScalar>();
-      E.block(2*F.rows(),1,F.rows(),1) = F.col(1).template cast<EScalar>();
-      return;
-  }
+  return oriented_facets(F,E);
 }
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
-template void igl::all_edges<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
-template void igl::all_edges<Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 2, 0, -1, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&);
-template void igl::all_edges<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 2, 0, -1, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&);
-template void igl::all_edges<Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
-template void igl::all_edges<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 2, 0, -1, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> >&);
+template void igl::all_edges<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::all_edges<Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 2, 0, -1, 2> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&);
+template void igl::all_edges<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 2, 0, -1, 2> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&);
+template void igl::all_edges<Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::all_edges<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 2, 0, -1, 2> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> >&);
 #endif

include/igl/all_edges.h

@@ -11,7 +11,10 @@
 #include <Eigen/Dense>
 namespace igl
 {
-  // ALL_EDGES Determines all "directed edges" of a given set of simplices
+  // Deprecated: call oriented_facets instead.
+  //
+  // ALL_EDGES Determines all "directed edges" of a given set of simplices. For
+  // a manifold mesh, this computes all of the half-edges
   //
   // Inputs:
   //   F  #F by simplex_size list of "faces"
@@ -24,7 +27,7 @@ namespace igl
   // once for each direction).
   template <typename DerivedF, typename DerivedE>
   IGL_INLINE void all_edges(
-    const Eigen::PlainObjectBase<DerivedF> & F,
+    const Eigen::MatrixBase<DerivedF> & F,
     Eigen::PlainObjectBase<DerivedE> & E);
 }
 

include/igl/biharmonic_coordinates.cpp

@@ -7,8 +7,11 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "biharmonic_coordinates.h"
 #include "cotmatrix.h"
+#include "sum.h"
 #include "massmatrix.h"
 #include "min_quad_with_fixed.h"
+#include "crouzeix_raviart_massmatrix.h"
+#include "crouzeix_raviart_cotmatrix.h"
 #include "normal_derivative.h"
 #include "on_boundary.h"
 #include <Eigen/Sparse>
@@ -45,11 +48,18 @@ IGL_INLINE bool igl::biharmonic_coordinates(
   // Subspace Design for Real-Time Shape Deformation" [Wang et al. 2015].
   SparseMatrix<double> A;
   {
-    SparseMatrix<double> N,Z,L,K,M;
-    normal_derivative(V,T,N);
-    Array<bool,Dynamic,1> I;
+    DiagonalMatrix<double,Dynamic> Minv;
+    SparseMatrix<double> L,K;
     Array<bool,Dynamic,Dynamic> C;
-    on_boundary(T,I,C);
+    {
+      Array<bool,Dynamic,1> I;
+      on_boundary(T,I,C);
+    }
+#ifdef false 
+    // Version described in paper is "wrong"
+    // http://www.cs.toronto.edu/~jacobson/images/error-in-linear-subspace-design-for-real-time-shape-deformation-2017-wang-et-al.pdf
+    SparseMatrix<double> N,Z,M;
+    normal_derivative(V,T,N);
     {
       std::vector<Triplet<double> >ZIJV;
       for(int t =0;t<T.rows();t++)
@@ -75,8 +85,50 @@ IGL_INLINE bool igl::biharmonic_coordinates(
     massmatrix(V,T,MASSMATRIX_TYPE_DEFAULT,M);
     // normalize
     M /= ((VectorXd)M.diagonal()).array().abs().maxCoeff();
-    DiagonalMatrix<double,Dynamic> Minv =
+    Minv =
       ((VectorXd)M.diagonal().array().inverse()).asDiagonal();
+#else
+    Eigen::SparseMatrix<double> M;
+    Eigen::MatrixXi E;
+    Eigen::VectorXi EMAP;
+    crouzeix_raviart_massmatrix(V,T,M,E,EMAP);
+    crouzeix_raviart_cotmatrix(V,T,E,EMAP,L);
+    // Ad  #E by #V facet-vertex incidence matrix
+    Eigen::SparseMatrix<double> Ad(E.rows(),V.rows());
+    {
+      std::vector<Eigen::Triplet<double> > AIJV(E.size());
+      for(int e = 0;e<E.rows();e++)
+      {
+        for(int c = 0;c<E.cols();c++)
+        {
+          AIJV[e+c*E.rows()] = Eigen::Triplet<double>(e,E(e,c),1);
+        }
+      }
+      Ad.setFromTriplets(AIJV.begin(),AIJV.end());
+    }
+    // Degrees
+    Eigen::VectorXd De;
+    sum(Ad,2,De);
+    Eigen::DiagonalMatrix<double,Eigen::Dynamic> De_diag = 
+      De.array().inverse().matrix().asDiagonal();
+    K = L*(De_diag*Ad);
+    // normalize
+    M /= ((VectorXd)M.diagonal()).array().abs().maxCoeff();
+    Minv = ((VectorXd)M.diagonal().array().inverse()).asDiagonal();
+    // kill boundary edges
+    for(int f = 0;f<T.rows();f++)
+    {
+      for(int c = 0;c<T.cols();c++)
+      {
+        if(C(f,c))
+        {
+          const int e = EMAP(f+T.rows()*c);
+          Minv.diagonal()(e) = 0;
+        }
+      }
+    }
+
+#endif
     switch(k)
     {
       default:

include/igl/copyleft/cgal/CSGTree.h

@@ -150,9 +150,9 @@ namespace igl
           // Returns mesh vertices in the desired output type, casting when
           // appropriate to floating precision.
           template <typename DerivedV>
-          Eigen::PlainObjectBase<DerivedV> cast_V() const
+          DerivedV cast_V() const
           {
-            Eigen::PlainObjectBase<DerivedV> dV;
+            DerivedV dV;
             dV.resize(m_V.rows(),m_V.cols());
             for(int i = 0;i<m_V.rows();i++)
             {

include/igl/copyleft/cgal/assign.cpp

@@ -10,7 +10,7 @@
 
 template <typename DerivedC, typename DerivedD>
 IGL_INLINE void igl::copyleft::cgal::assign(
-  const Eigen::PlainObjectBase<DerivedC> & C,
+  const Eigen::MatrixBase<DerivedC> & C,
   Eigen::PlainObjectBase<DerivedD> & D)
 {
   D.resizeLike(C);
@@ -23,32 +23,42 @@ IGL_INLINE void igl::copyleft::cgal::assign(
   }
 }
 
+template <typename ReturnScalar, typename DerivedC>
+IGL_INLINE 
+Eigen::Matrix<
+  ReturnScalar,
+  DerivedC::RowsAtCompileTime, 
+  DerivedC::ColsAtCompileTime, 
+  1,
+  DerivedC::MaxRowsAtCompileTime, 
+  DerivedC::MaxColsAtCompileTime>
+igl::copyleft::cgal::assign(
+  const Eigen::MatrixBase<DerivedC> & C)
+{
+  Eigen::Matrix<
+    ReturnScalar,
+    DerivedC::RowsAtCompileTime, 
+    DerivedC::ColsAtCompileTime, 
+    1,
+    DerivedC::MaxRowsAtCompileTime, 
+    DerivedC::MaxColsAtCompileTime> D;
+  assign(C,D);
+  return D;
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
-// generated by autoexplicit.sh
-template void igl::copyleft::cgal::assign<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
-// generated by autoexplicit.sh
-template void igl::copyleft::cgal::assign<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&);
-// generated by autoexplicit.sh
-template void igl::copyleft::cgal::assign<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
-// generated by autoexplicit.sh
-template void igl::copyleft::cgal::assign<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> >&);
-// generated by autoexplicit.sh
-template void igl::copyleft::cgal::assign<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
-// generated by autoexplicit.sh
-template void igl::copyleft::cgal::assign<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> >&);
-// generated by autoexplicit.sh
-template void igl::copyleft::cgal::assign<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
-// generated by autoexplicit.sh
-template void igl::copyleft::cgal::assign<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&);
-// generated by autoexplicit.sh
-template void igl::copyleft::cgal::assign<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
-// generated by autoexplicit.sh
-template void igl::copyleft::cgal::assign<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&);
-// generated by autoexplicit.sh
-template void igl::copyleft::cgal::assign<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
-// generated by autoexplicit.sh
-template void igl::copyleft::cgal::assign<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, 8, 3, 0, 8, 3>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, 8, 3, 0, 8, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, 8, 3, 0, 8, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, 8, 3, 0, 8, 3> >&);
-// generated by autoexplicit.sh
-template void igl::copyleft::cgal::assign<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, 8, 3, 0, 8, 3>, Eigen::Matrix<double, 8, 3, 0, 8, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, 8, 3, 0, 8, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 8, 3, 0, 8, 3> >&);
+template void igl::copyleft::cgal::assign<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
+template void igl::copyleft::cgal::assign<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&);
+template void igl::copyleft::cgal::assign<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::copyleft::cgal::assign<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> >&);
+template void igl::copyleft::cgal::assign<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::copyleft::cgal::assign<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> >&);
+template void igl::copyleft::cgal::assign<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
+template void igl::copyleft::cgal::assign<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&);
+template void igl::copyleft::cgal::assign<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::copyleft::cgal::assign<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&);
+template void igl::copyleft::cgal::assign<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
+template void igl::copyleft::cgal::assign<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, 8, 3, 0, 8, 3>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, 8, 3, 0, 8, 3> >(Eigen::MatrixBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, 8, 3, 0, 8, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, 8, 3, 0, 8, 3> >&);
+template void igl::copyleft::cgal::assign<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, 8, 3, 0, 8, 3>, Eigen::Matrix<double, 8, 3, 0, 8, 3> >(Eigen::MatrixBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, 8, 3, 0, 8, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 8, 3, 0, 8, 3> >&);
 #endif

include/igl/copyleft/cgal/assign.h

@@ -19,8 +19,19 @@ namespace igl
     {
       template <typename DerivedC, typename DerivedD>
       IGL_INLINE void assign(
-        const Eigen::PlainObjectBase<DerivedC> & C,
+        const Eigen::MatrixBase<DerivedC> & C,
         Eigen::PlainObjectBase<DerivedD> & D);
+      template <typename ReturnScalar, typename DerivedC>
+      IGL_INLINE 
+      Eigen::Matrix<
+        ReturnScalar,
+        DerivedC::RowsAtCompileTime, 
+        DerivedC::ColsAtCompileTime, 
+        1,
+        DerivedC::MaxRowsAtCompileTime, 
+        DerivedC::MaxColsAtCompileTime> 
+      assign(
+        const Eigen::MatrixBase<DerivedC> & C);
     }
   }
 }

include/igl/copyleft/cgal/assign_scalar.cpp

@@ -64,3 +64,36 @@ IGL_INLINE void igl::copyleft::cgal::assign_scalar(
 {
   d = c;
 }
+
+IGL_INLINE void igl::copyleft::cgal::assign_scalar(
+  const CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt::FT & cgal,
+  CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt::FT & d)
+{
+  d = cgal;
+}
+
+IGL_INLINE void igl::copyleft::cgal::assign_scalar(
+  const CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt::FT & cgal,
+  double & d)
+{
+  const auto interval = CGAL::to_interval(cgal);
+  d = interval.first;
+  do {
+      const double next = nextafter(d, interval.second);
+      if (CGAL::abs(cgal-d) < CGAL::abs(cgal-next)) break;
+      d = next;
+  } while (d < interval.second);
+}
+
+IGL_INLINE void igl::copyleft::cgal::assign_scalar(
+  const CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt::FT & cgal,
+  float& d)
+{
+  const auto interval = CGAL::to_interval(cgal);
+  d = interval.first;
+  do {
+      const float next = nextafter(d, float(interval.second));
+      if (CGAL::abs(cgal-d) < CGAL::abs(cgal-next)) break;
+      d = next;
+  } while (d < float(interval.second));
+}

include/igl/copyleft/cgal/assign_scalar.h

@@ -8,7 +8,8 @@
 #ifndef IGL_COPYLEFT_CGAL_ASSIGN_SCALAR_H
 #define IGL_COPYLEFT_CGAL_ASSIGN_SCALAR_H
 #include "../../igl_inline.h"
-#include "CGAL_includes.hpp"
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel_with_sqrt.h>
 namespace igl
 {
   namespace copyleft
@@ -37,6 +38,15 @@ namespace igl
       IGL_INLINE void assign_scalar(
         const float& c,
         double& d);
+      IGL_INLINE void assign_scalar(
+        const CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt::FT & cgal,
+        CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt::FT & d);
+      IGL_INLINE void assign_scalar(
+        const CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt::FT & cgal,
+        double & d);
+      IGL_INLINE void assign_scalar(
+        const CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt::FT & cgal,
+        float& d);
     }
   }
 }

include/igl/copyleft/cgal/remesh_intersections.cpp

@@ -448,7 +448,9 @@ IGL_INLINE void igl::copyleft::cgal::remesh_intersections(
       std::transform(FF.data(), FF.data() + FF.rows()*FF.cols(),
           FF.data(), [&vv_to_unique](const typename DerivedFF::Scalar& a)
           { return vv_to_unique[a]; });
-      IM = igl::LinSpaced<
+      IM.resize(unique_vv.rows());
+      // Have to use << instead of = becasue Eigen's PlainObjectBase is annoying
+      IM << igl::LinSpaced<
         Eigen::Matrix<typename DerivedIM::Scalar, Eigen::Dynamic,1 >
         >(unique_vv.rows(), 0, unique_vv.rows()-1);
     }else 

include/igl/copyleft/tetgen/tetgenio_to_tetmesh.cpp

@@ -136,7 +136,7 @@ IGL_INLINE bool igl::copyleft::tetgen::tetgenio_to_tetmesh(
   Eigen::PlainObjectBase<DerivedV>& V,
   Eigen::PlainObjectBase<DerivedT>& T)
 {
-  Eigen::PlainObjectBase<DerivedT> F;
+  Eigen::Matrix<typename DerivedT::Scalar,Eigen::Dynamic,3> F;
   return tetgenio_to_tetmesh(out,V,T,F);
 }
 

include/igl/crouzeix_raviart_cotmatrix.cpp

@@ -0,0 +1,100 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2017 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 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "crouzeix_raviart_cotmatrix.h"
+#include "unique_simplices.h"
+#include "oriented_facets.h"
+#include "is_edge_manifold.h"
+#include "cotmatrix_entries.h"
+
+template <typename DerivedV, typename DerivedF, typename LT, typename DerivedE, typename DerivedEMAP>
+void igl::crouzeix_raviart_cotmatrix(
+  const Eigen::MatrixBase<DerivedV> & V, 
+  const Eigen::MatrixBase<DerivedF> & F, 
+  Eigen::SparseMatrix<LT> & L,
+  Eigen::PlainObjectBase<DerivedE> & E,
+  Eigen::PlainObjectBase<DerivedEMAP> & EMAP)
+{
+  // All occurances of directed "facets"
+  Eigen::MatrixXi allE;
+  oriented_facets(F,allE);
+  Eigen::VectorXi _1;
+  unique_simplices(allE,E,_1,EMAP);
+  return crouzeix_raviart_cotmatrix(V,F,E,EMAP,L);
+}
+
+template <typename DerivedV, typename DerivedF, typename DerivedE, typename DerivedEMAP, typename LT>
+void igl::crouzeix_raviart_cotmatrix(
+  const Eigen::MatrixBase<DerivedV> & V, 
+  const Eigen::MatrixBase<DerivedF> & F, 
+  const Eigen::MatrixBase<DerivedE> & E,
+  const Eigen::MatrixBase<DerivedEMAP> & EMAP,
+  Eigen::SparseMatrix<LT> & L)
+{
+  // number of rows
+  const int m = F.rows();
+  // Element simplex size
+  const int ss = F.cols();
+  // Mesh should be edge-manifold
+  assert(F.cols() != 3 || is_edge_manifold(F));
+  typedef Eigen::Matrix<LT,Eigen::Dynamic,Eigen::Dynamic> MatrixXS;
+  MatrixXS C;
+  cotmatrix_entries(V,F,C);
+  Eigen::MatrixXi F2E(m,ss);
+  {
+    int k =0;
+    for(int c = 0;c<ss;c++)
+    {
+      for(int f = 0;f<m;f++)
+      {
+        F2E(f,c) = k++;
+      }
+    }
+  }
+  // number of entries inserted per facet
+  const int k = ss*(ss-1)*2;
+  std::vector<Eigen::Triplet<LT> > LIJV;LIJV.reserve(k*m);
+  Eigen::VectorXi LI(k),LJ(k),LV(k);
+  // Compensation factor to match scales in matlab version
+  double factor = 2.0;
+
+  switch(ss)
+  {
+    default: assert(false && "unsupported simplex size");
+    case 3:
+      factor = 4.0;
+      LI<<0,1,2,1,2,0,0,1,2,1,2,0;
+      LJ<<1,2,0,0,1,2,0,1,2,1,2,0;
+      LV<<2,0,1,2,0,1,2,0,1,2,0,1;
+      break;
+    case 4:
+      factor *= -1.0;
+      LI<<0,3,3,3,1,2,1,0,1,2,2,0,0,3,3,3,1,2,1,0,1,2,2,0;
+      LJ<<1,0,1,2,2,0,0,3,3,3,1,2,0,3,3,3,1,2,1,0,1,2,2,0;
+      LV<<2,3,4,5,0,1,2,3,4,5,0,1,2,3,4,5,0,1,2,3,4,5,0,1;
+      break;
+  }
+
+  for(int f=0;f<m;f++)
+  {
+    for(int c = 0;c<k;c++)
+    {
+      LIJV.emplace_back(
+        EMAP(F2E(f,LI(c))),
+        EMAP(F2E(f,LJ(c))),
+        (c<(k/2)?-1.:1.) * factor *C(f,LV(c)));
+    }
+  }
+  L.resize(E.rows(),E.rows());
+  L.setFromTriplets(LIJV.begin(),LIJV.end());
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::crouzeix_raviart_cotmatrix<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::SparseMatrix<double, 0, int>&);
+#endif

include/igl/crouzeix_raviart_cotmatrix.h

@@ -0,0 +1,49 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2017 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 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_CROUZEIX_RAVIART_COTMATRIX
+#define IGL_CROUZEIX_RAVIART_COTMATRIX
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <Eigen/Sparse>
+namespace igl
+{
+  // CROUZEIX_RAVIART_COTMATRIX Compute the Crouzeix-Raviart cotangent
+  // stiffness matrix.
+  //
+  // See for example "Discrete Quadratic Curvature Energies" [Wardetzky, Bergou,
+  // Harmon, Zorin, Grinspun 2007]
+  //
+  // Inputs:
+  //   V  #V by dim list of vertex positions
+  //   F  #F by 3/4 list of triangle/tetrahedron indices
+  // Outputs:
+  //   L  #E by #E edge/face-based diagonal cotangent matrix
+  //   E  #E by 2/3 list of edges/faces
+  //   EMAP  #F*3/4 list of indices mapping allE to E
+  //
+  // See also: crouzeix_raviart_massmatrix
+  template <typename DerivedV, typename DerivedF, typename LT, typename DerivedE, typename DerivedEMAP>
+  void crouzeix_raviart_cotmatrix(
+      const Eigen::MatrixBase<DerivedV> & V, 
+      const Eigen::MatrixBase<DerivedF> & F, 
+      Eigen::SparseMatrix<LT> & L,
+      Eigen::PlainObjectBase<DerivedE> & E,
+      Eigen::PlainObjectBase<DerivedEMAP> & EMAP);
+  // wrapper if E and EMAP are already computed (better match!)
+  template <typename DerivedV, typename DerivedF, typename DerivedE, typename DerivedEMAP, typename LT>
+  void crouzeix_raviart_cotmatrix(
+      const Eigen::MatrixBase<DerivedV> & V, 
+      const Eigen::MatrixBase<DerivedF> & F, 
+      const Eigen::MatrixBase<DerivedE> & E,
+      const Eigen::MatrixBase<DerivedEMAP> & EMAP,
+      Eigen::SparseMatrix<LT> & L);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "crouzeix_raviart_cotmatrix.cpp"
+#endif
+#endif

include/igl/crouzeix_raviart_massmatrix.cpp

@@ -7,25 +7,26 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "crouzeix_raviart_massmatrix.h"
 #include "unique_simplices.h"
-#include "all_edges.h"
+#include "oriented_facets.h"
 
 #include "is_edge_manifold.h"
 #include "doublearea.h"
+#include "volume.h"
 
 #include <cassert>
 #include <vector>
 
 template <typename MT, typename DerivedV, typename DerivedF, typename DerivedE, typename DerivedEMAP>
 void igl::crouzeix_raviart_massmatrix(
-    const Eigen::PlainObjectBase<DerivedV> & V, 
-    const Eigen::PlainObjectBase<DerivedF> & F, 
+    const Eigen::MatrixBase<DerivedV> & V, 
+    const Eigen::MatrixBase<DerivedF> & F, 
     Eigen::SparseMatrix<MT> & M,
     Eigen::PlainObjectBase<DerivedE> & E,
     Eigen::PlainObjectBase<DerivedEMAP> & EMAP)
 {
-  // All occurances of directed edges
+  // All occurances of directed "facets"
   Eigen::MatrixXi allE;
-  all_edges(F,allE);
+  oriented_facets(F,allE);
   Eigen::VectorXi _1;
   unique_simplices(allE,E,_1,EMAP);
   return crouzeix_raviart_massmatrix(V,F,E,EMAP,M);
@@ -33,28 +34,42 @@ void igl::crouzeix_raviart_massmatrix(
 
 template <typename MT, typename DerivedV, typename DerivedF, typename DerivedE, typename DerivedEMAP>
 void igl::crouzeix_raviart_massmatrix(
-    const Eigen::PlainObjectBase<DerivedV> & V, 
-    const Eigen::PlainObjectBase<DerivedF> & F, 
-    const Eigen::PlainObjectBase<DerivedE> & E,
-    const Eigen::PlainObjectBase<DerivedEMAP> & EMAP,
+    const Eigen::MatrixBase<DerivedV> & V, 
+    const Eigen::MatrixBase<DerivedF> & F, 
+    const Eigen::MatrixBase<DerivedE> & E,
+    const Eigen::MatrixBase<DerivedEMAP> & EMAP,
     Eigen::SparseMatrix<MT> & M)
 {
   using namespace Eigen;
   using namespace std;
-  assert(F.cols() == 3);
-  // Mesh should be edge-manifold
-  assert(is_edge_manifold(F));
+  // Mesh should be edge-manifold (TODO: replace `is_edge_manifold` with
+  // `is_facet_manifold`)
+  assert(F.cols() != 3 || is_edge_manifold(F));
   // number of elements (triangles)
-  int m = F.rows();
-  // Get triangle areas
+  const int m = F.rows();
+  // Get triangle areas/volumes
   VectorXd TA;
-  doublearea(V,F,TA);
-  vector<Triplet<MT> > MIJV(3*m);
+  // Element simplex size
+  const int ss = F.cols();
+  switch(ss)
+  {
+    default:
+      assert(false && "Unsupported simplex size");
+    case 3:
+      doublearea(V,F,TA);
+      TA *= 0.5;
+      break;
+    case 4:
+      volume(V,F,TA);
+      break;
+  }
+  vector<Triplet<MT> > MIJV(ss*m);
+  assert(EMAP.size() == m*ss);
   for(int f = 0;f<m;f++)
   {
-    for(int c = 0;c<3;c++)
+    for(int c = 0;c<ss;c++)
     {
-      MIJV[f+m*c] = Triplet<MT>(EMAP(f+m*c),EMAP(f+m*c),TA(f)*0.5);
+      MIJV[f+m*c] = Triplet<MT>(EMAP(f+m*c),EMAP(f+m*c),TA(f)/(double)(ss));
     }
   }
   M.resize(E.rows(),E.rows());
@@ -63,6 +78,8 @@ void igl::crouzeix_raviart_massmatrix(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
-template void igl::crouzeix_raviart_massmatrix<double, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::SparseMatrix<double, 0, int>&);
-template void igl::crouzeix_raviart_massmatrix<float, Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::SparseMatrix<float, 0, int>&);
+// generated by autoexplicit.sh
+template void igl::crouzeix_raviart_massmatrix<double, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int>&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template void igl::crouzeix_raviart_massmatrix<double, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::SparseMatrix<double, 0, int>&);
+template void igl::crouzeix_raviart_massmatrix<float, Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::SparseMatrix<float, 0, int>&);
 #endif

include/igl/crouzeix_raviart_massmatrix.h

@@ -19,38 +19,29 @@ namespace igl
   // See for example "Discrete Quadratic Curvature Energies" [Wardetzky, Bergou,
   // Harmon, Zorin, Grinspun 2007]
   //
-  // Templates:
-  //   MT  type of eigen sparse matrix for M (e.g. double for
-  //     SparseMatrix<double>)
-  //   DerivedV  derived type of eigen matrix for V (e.g. derived from
-  //     MatrixXd)
-  //   DerivedF  derived type of eigen matrix for F (e.g. derived from
-  //     MatrixXi)
-  //   DerivedE  derived type of eigen matrix for E (e.g. derived from
-  //     MatrixXi)
   // Inputs:
   //   V  #V by dim list of vertex positions
-  //   F  #F by 3 list of triangle indices
+  //   F  #F by 3/4 list of triangle/tetrahedron indices
   // Outputs:
-  //   M  #E by #E edge-based diagonal mass matrix
-  //   E  #E by 2 list of edges
-  //   EMAP  #F*3 list of indices mapping allE to E
-  //
+  //   M  #E by #E edge/face-based diagonal mass matrix
+  //   E  #E by 2/3 list of edges/faces
+  //   EMAP  #F*3/4 list of indices mapping allE to E
   //
+  // See also: crouzeix_raviart_cotmatrix
   template <typename MT, typename DerivedV, typename DerivedF, typename DerivedE, typename DerivedEMAP>
   void crouzeix_raviart_massmatrix(
-      const Eigen::PlainObjectBase<DerivedV> & V, 
-      const Eigen::PlainObjectBase<DerivedF> & F, 
+      const Eigen::MatrixBase<DerivedV> & V, 
+      const Eigen::MatrixBase<DerivedF> & F, 
       Eigen::SparseMatrix<MT> & M,
       Eigen::PlainObjectBase<DerivedE> & E,
       Eigen::PlainObjectBase<DerivedEMAP> & EMAP);
   // wrapper if E and EMAP are already computed (better match!)
   template <typename MT, typename DerivedV, typename DerivedF, typename DerivedE, typename DerivedEMAP>
   void crouzeix_raviart_massmatrix(
-      const Eigen::PlainObjectBase<DerivedV> & V, 
-      const Eigen::PlainObjectBase<DerivedF> & F, 
-      const Eigen::PlainObjectBase<DerivedE> & E,
-      const Eigen::PlainObjectBase<DerivedEMAP> & EMAP,
+      const Eigen::MatrixBase<DerivedV> & V, 
+      const Eigen::MatrixBase<DerivedF> & F, 
+      const Eigen::MatrixBase<DerivedE> & E,
+      const Eigen::MatrixBase<DerivedEMAP> & EMAP,
       Eigen::SparseMatrix<MT> & M);
 }
 #ifndef IGL_STATIC_LIBRARY

include/igl/cut_mesh.cpp

@@ -12,6 +12,7 @@
 #include <igl/HalfEdgeIterator.h>
 #include <set>
 
+// This file violates many of the libigl style guidelines.
 
 namespace igl {
 
@@ -22,15 +23,16 @@ namespace igl {
   public:
     // Input
     //mesh
-    const Eigen::PlainObjectBase<DerivedV> &V;
-    const Eigen::PlainObjectBase<DerivedF> &F;
-    const Eigen::PlainObjectBase<DerivedTT> &TT;
-    const Eigen::PlainObjectBase<DerivedTT> &TTi;
+    const DerivedV &V;
+    const DerivedF &F;
+    // TT is the same type as TTi? This is likely to break at some point
+    const DerivedTT &TT;
+    const DerivedTT &TTi;
     const std::vector<std::vector<VFType> >& VF;
     const std::vector<std::vector<VFType> >& VFi;
     const std::vector<bool> &V_border; // bool
     //edges to cut
-    const Eigen::PlainObjectBase<DerivedC> &Handle_Seams; // 3 bool
+    const DerivedC &Handle_Seams; // 3 bool
 
     // total number of scalar variables
     int num_scalar_variables;
@@ -41,14 +43,15 @@ namespace igl {
     // per vertex variable indexes
     std::vector<std::vector<int> > HandleV_Integer;
 
-    IGL_INLINE MeshCutterMini(const Eigen::PlainObjectBase<DerivedV> &_V,
-                              const Eigen::PlainObjectBase<DerivedF> &_F,
-                              const Eigen::PlainObjectBase<DerivedTT> &_TT,
-                              const Eigen::PlainObjectBase<DerivedTT> &_TTi,
-                              const std::vector<std::vector<VFType> > &_VF,
-                              const std::vector<std::vector<VFType> > &_VFi,
-                              const std::vector<bool> &_V_border,
-                              const Eigen::PlainObjectBase<DerivedC> &_Handle_Seams);
+    IGL_INLINE MeshCutterMini(
+      const Eigen::PlainObjectBase<DerivedV> &_V,
+      const Eigen::PlainObjectBase<DerivedF> &_F,
+      const Eigen::PlainObjectBase<DerivedTT> &_TT,
+      const Eigen::PlainObjectBase<DerivedTT> &_TTi,
+      const std::vector<std::vector<VFType> > &_VF,
+      const std::vector<std::vector<VFType> > &_VFi,
+      const std::vector<bool> &_V_border,
+      const Eigen::PlainObjectBase<DerivedC> &_Handle_Seams);
 
     // vertex to variable mapping
     // initialize the mapping for a given sampled mesh
@@ -80,27 +83,26 @@ namespace igl {
 
 template <typename DerivedV, typename DerivedF, typename VFType, typename DerivedTT, typename DerivedC>
 IGL_INLINE igl::MeshCutterMini<DerivedV, DerivedF, VFType, DerivedTT, DerivedC>::
-MeshCutterMini(const Eigen::PlainObjectBase<DerivedV> &_V,
-               const Eigen::PlainObjectBase<DerivedF> &_F,
-               const Eigen::PlainObjectBase<DerivedTT> &_TT,
-               const Eigen::PlainObjectBase<DerivedTT> &_TTi,
-               const std::vector<std::vector<VFType> > &_VF,
-               const std::vector<std::vector<VFType> > &_VFi,
-               const std::vector<bool> &_V_border,
-               const Eigen::PlainObjectBase<DerivedC> &_Handle_Seams):
-V(_V),
-F(_F),
-TT(_TT),
-TTi(_TTi),
-VF(_VF),
-VFi(_VFi),
-V_border(_V_border),
-Handle_Seams(_Handle_Seams)
+MeshCutterMini(
+  const Eigen::PlainObjectBase<DerivedV> &_V,
+  const Eigen::PlainObjectBase<DerivedF> &_F,
+  const Eigen::PlainObjectBase<DerivedTT> &_TT,
+  const Eigen::PlainObjectBase<DerivedTT> &_TTi,
+  const std::vector<std::vector<VFType> > &_VF,
+  const std::vector<std::vector<VFType> > &_VFi,
+  const std::vector<bool> &_V_border,
+  const Eigen::PlainObjectBase<DerivedC> &_Handle_Seams):
+  V(_V),
+  F(_F),
+  TT(_TT),
+  TTi(_TTi),
+  VF(_VF),
+  VFi(_VFi),
+  V_border(_V_border),
+  Handle_Seams(_Handle_Seams)
 {
   num_scalar_variables=0;
-
   HandleS_Index.setConstant(F.rows(),3,-1);
-
   HandleV_Integer.resize(V.rows());
 }
 
@@ -151,7 +153,7 @@ FindInitialPos(const int vert,
   int f_init;
   int edge_init;
   FirstPos(vert,f_init,edge_init); // todo manually the function
-  igl::HalfEdgeIterator<DerivedF> VFI(F,TT,TTi,f_init,edge_init);
+  igl::HalfEdgeIterator<DerivedF,DerivedTT,DerivedTT> VFI(F,TT,TTi,f_init,edge_init);
 
   bool vertexB = V_border[vert];
   bool possible_split=false;
@@ -198,7 +200,7 @@ MapIndexes(const int  vert,
   ///insert an initial index
   int curr_index=AddNewIndex(vert);
   ///and initialize the jumping pos
-  igl::HalfEdgeIterator<DerivedF> VFI(F,TT,TTi,f_init,edge_init);
+  igl::HalfEdgeIterator<DerivedF,DerivedTT,DerivedTT> VFI(F,TT,TTi,f_init,edge_init);
   bool complete_turn=false;
   do
   {
@@ -235,7 +237,7 @@ InitMappingSeam(const int vert)
   int f_init = VF[vert][0];
   int indexE = VFi[vert][0];
 
-  igl::HalfEdgeIterator<DerivedF> VFI(F,TT,TTi,f_init,indexE);
+  igl::HalfEdgeIterator<DerivedF,DerivedTT,DerivedTT> VFI(F,TT,TTi,f_init,indexE);
 
   int edge_init;
   int face_init;
@@ -260,16 +262,16 @@ InitMappingSeam()
 
 template <typename DerivedV, typename DerivedF, typename VFType, typename DerivedTT, typename DerivedC>
 IGL_INLINE void igl::cut_mesh(
-                                                                  const Eigen::PlainObjectBase<DerivedV> &V,
-                                                                  const Eigen::PlainObjectBase<DerivedF> &F,
-                                                                  const std::vector<std::vector<VFType> >& VF,
-                                                                  const std::vector<std::vector<VFType> >& VFi,
-                                                                  const Eigen::PlainObjectBase<DerivedTT>& TT,
-                                                                  const Eigen::PlainObjectBase<DerivedTT>& TTi,
-                                                                  const std::vector<bool> &V_border,
-                                                                  const Eigen::PlainObjectBase<DerivedC> &cuts,
-                                                                  Eigen::PlainObjectBase<DerivedV> &Vcut,
-                                                                  Eigen::PlainObjectBase<DerivedF> &Fcut)
+  const Eigen::PlainObjectBase<DerivedV> &V,
+  const Eigen::PlainObjectBase<DerivedF> &F,
+  const std::vector<std::vector<VFType> >& VF,
+  const std::vector<std::vector<VFType> >& VFi,
+  const Eigen::PlainObjectBase<DerivedTT>& TT,
+  const Eigen::PlainObjectBase<DerivedTT>& TTi,
+  const std::vector<bool> &V_border,
+  const Eigen::PlainObjectBase<DerivedC> &cuts,
+  Eigen::PlainObjectBase<DerivedV> &Vcut,
+  Eigen::PlainObjectBase<DerivedF> &Fcut)
 {
   //finding the cuts is done, now we need to actually generate a cut mesh
   igl::MeshCutterMini<DerivedV, DerivedF, VFType, DerivedTT, DerivedC> mc(V, F, TT, TTi, VF, VFi, V_border, cuts);
@@ -302,26 +304,21 @@ IGL_INLINE void igl::cut_mesh(
 //Wrapper of the above with only vertices and faces as mesh input
 template <typename DerivedV, typename DerivedF, typename DerivedC>
 IGL_INLINE void igl::cut_mesh(
-                                                                  const Eigen::PlainObjectBase<DerivedV> &V,
-                                                                  const Eigen::PlainObjectBase<DerivedF> &F,
-                                                                  const Eigen::PlainObjectBase<DerivedC> &cuts,
-                                                                  Eigen::PlainObjectBase<DerivedV> &Vcut,
-                                                                  Eigen::PlainObjectBase<DerivedF> &Fcut)
+  const Eigen::PlainObjectBase<DerivedV> &V,
+  const Eigen::PlainObjectBase<DerivedF> &F,
+  const Eigen::PlainObjectBase<DerivedC> &cuts,
+  Eigen::PlainObjectBase<DerivedV> &Vcut,
+  Eigen::PlainObjectBase<DerivedF> &Fcut)
 {
-
   std::vector<std::vector<int> > VF, VFi;
   igl::vertex_triangle_adjacency(V,F,VF,VFi);
-
+  // Alec: Cast? Why? This is likely to break.
   Eigen::MatrixXd Vt = V;
   Eigen::MatrixXi Ft = F;
   Eigen::MatrixXi TT, TTi;
   igl::triangle_triangle_adjacency(Ft,TT,TTi);
-
   std::vector<bool> V_border = igl::is_border_vertex(V,F);
-
   igl::cut_mesh(V, F, VF, VFi, TT, TTi, V_border, cuts, Vcut, Fcut);
-
-
 }
 
 #ifdef IGL_STATIC_LIBRARY

include/igl/cut_mesh.h

@@ -20,6 +20,8 @@ namespace igl
   // Todo: this combinatorial operation should not depend on the vertex
   // positions V.
   //
+  // Known issues: Assumes mesh is edge-manifold.
+  //
   // Inputs:
   //   V  #V by 3 list of the vertex positions
   //   F  #F by 3 list of the faces (must be triangles)

include/igl/eigs.cpp

@@ -45,8 +45,10 @@ IGL_INLINE bool igl::eigs(
   Scalar conv = 1e-14;
   int max_iter = 100;
   int i = 0;
+  //std::cout<<"start"<<std::endl;
   while(true)
   {
+    //std::cout<<i<<std::endl;
     // Random initial guess
     VectorXS y = VectorXS::Random(n,1);
     Scalar eff_sigma = 0;
@@ -132,8 +134,13 @@ IGL_INLINE bool igl::eigs(
       cerr<<"Failed to converge."<<endl;
       return false;
     }
-    if(i==0 || (S.head(i).array()-sigma).abs().maxCoeff()>1e-14)
+    if(
+      i==0 || 
+      (S.head(i).array()-sigma).abs().maxCoeff()>1e-14 ||
+      ((U.leftCols(i).transpose()*B*x).array().abs()<=1e-7).all()
+      )
     {
+      //cout<<"Found "<<i<<"th mode"<<endl;
       U.col(i) = x;
       S(i) = sigma;
       i++;
@@ -143,14 +150,18 @@ IGL_INLINE bool igl::eigs(
       }
     }else
     {
+      //std::cout<<"i: "<<i<<std::endl;
+      //std::cout<<"  "<<S.head(i).transpose()<<" << "<<sigma<<std::endl;
+      //std::cout<<"  "<<(S.head(i).array()-sigma).abs().maxCoeff()<<std::endl;
+      //std::cout<<"  "<<(U.leftCols(i).transpose()*B*x).array().abs().transpose()<<std::endl;
       // restart with new random guess.
-      cout<<"RESTART!"<<endl;
+      cout<<"igl::eigs RESTART"<<endl;
     }
   }
   // finally sort
   VectorXi I;
   igl::sort(S,1,false,sS,I);
-  sU = igl::slice(U,I,2);
+  igl::slice(U,I,2,sU);
   sS /= rescale;
   sU /= sqrt(rescale);
   return true;

include/igl/eigs.h

@@ -27,12 +27,12 @@ namespace igl
   //   A  #A by #A symmetric matrix
   //   B  #A by #A symmetric positive-definite matrix
   //   k  number of eigen pairs to compute
+  //   type  whether to extract from the high or low end
   // Outputs:
   //   sU  #A by k list of sorted eigen vectors (descending)
   //   sS  k list of sorted eigen values (descending)
   //
   // Known issues:
-  //   - only one pair per eigen value is found (no multiples)
   //   - only the 'sm' small magnitude eigen values are well supported
   //   
   enum EigsType

include/igl/embree/bone_visible.cpp

@@ -25,7 +25,7 @@ IGL_INLINE void igl::embree::bone_visible(
   Eigen::PlainObjectBase<Derivedflag>  & flag)
 {
   // "double sided lighting"
-  Eigen::PlainObjectBase<DerivedF> FF;
+  Eigen::Matrix<typename DerivedF::Scalar,Eigen::Dynamic,Eigen::Dynamic> FF;
   FF.resize(F.rows()*2,F.cols());
   FF << F, F.rowwise().reverse();
   // Initialize intersector

include/igl/exterior_edges.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 "exterior_edges.h"
-#include "all_edges.h"
+#include "oriented_facets.h"
 #include "sort.h"
 #include "unique.h"
 
@@ -51,7 +51,7 @@ IGL_INLINE void igl::exterior_edges(
   const size_t m = F.rows();
   MatrixXi all_E,sall_E,sort_order;
   // Sort each edge by index
-  all_edges(F,all_E);
+  oriented_facets(F,all_E);
   sort(all_E,2,true,sall_E,sort_order);
   // Find unique edges
   MatrixXi uE;

include/igl/find.cpp

@@ -121,19 +121,18 @@ IGL_INLINE void igl::find(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
-// generated by autoexplicit.sh
-template void igl::find<Eigen::CwiseBinaryOp<Eigen::internal::scalar_cmp_op<long, (Eigen::internal::ComparisonName)0>, Eigen::PartialReduxExpr<Eigen::Array<bool, -1, 3, 0, -1, 3>, Eigen::internal::member_count<long>, 1> const, Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<long>, Eigen::Array<long, -1, 1, 0, -1, 1> > const>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::DenseBase<Eigen::CwiseBinaryOp<Eigen::internal::scalar_cmp_op<long, (Eigen::internal::ComparisonName)0>, Eigen::PartialReduxExpr<Eigen::Array<bool, -1, 3, 0, -1, 3>, Eigen::internal::member_count<long>, 1> const, Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<long>, Eigen::Array<long, -1, 1, 0, -1, 1> > const> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-// generated by autoexplicit.sh
+
 template void igl::find<bool, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Array<bool, -1, 1, 0, -1, 1> >(Eigen::SparseMatrix<bool, 0, int> const&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::DenseBase<Eigen::Array<bool, -1, 1, 0, -1, 1> >&);
-// generated by autoexplicit.sh
 template void igl::find<Eigen::Array<bool, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::DenseBase<Eigen::Array<bool, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-// generated by autoexplicit.sh
-template void igl::find<Eigen::CwiseBinaryOp<Eigen::internal::scalar_cmp_op<int, (Eigen::internal::ComparisonName)0>, Eigen::Array<int, -1, 1, 0, -1, 1> const, Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<int>, Eigen::Array<int, -1, 1, 0, -1, 1> > const>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::DenseBase<Eigen::CwiseBinaryOp<Eigen::internal::scalar_cmp_op<int, (Eigen::internal::ComparisonName)0>, Eigen::Array<int, -1, 1, 0, -1, 1> const, Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<int>, Eigen::Array<int, -1, 1, 0, -1, 1> > const> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-// generated by autoexplicit.sh
 template void igl::find<double, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::SparseMatrix<double, 0, int> const&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::DenseBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 template void igl::find<double, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::SparseMatrix<double, 0, int> const&, Eigen::DenseBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, Eigen::DenseBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, Eigen::DenseBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
-// generated by autoexplicit.sh
 template void igl::find<double, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::SparseMatrix<double, 0, int> const&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::DenseBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template void igl::find<Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::find<double, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::SparseMatrix<double, 0, int> const&, Eigen::DenseBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::DenseBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::DenseBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+#if EIGEN_VERSION_AT_LEAST(3,3,0)
+#else
+template void igl::find<Eigen::CwiseBinaryOp<Eigen::internal::scalar_cmp_op<long, (Eigen::internal::ComparisonName)0>, Eigen::PartialReduxExpr<Eigen::Array<bool, -1, 3, 0, -1, 3>, Eigen::internal::member_count<long>, 1> const, Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<long>, Eigen::Array<long, -1, 1, 0, -1, 1> > const>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::DenseBase<Eigen::CwiseBinaryOp<Eigen::internal::scalar_cmp_op<long, (Eigen::internal::ComparisonName)0>, Eigen::PartialReduxExpr<Eigen::Array<bool, -1, 3, 0, -1, 3>, Eigen::internal::member_count<long>, 1> const, Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<long>, Eigen::Array<long, -1, 1, 0, -1, 1> > const> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template void igl::find<Eigen::CwiseBinaryOp<Eigen::internal::scalar_cmp_op<int, (Eigen::internal::ComparisonName)0>, Eigen::Array<int, -1, 1, 0, -1, 1> const, Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<int>, Eigen::Array<int, -1, 1, 0, -1, 1> > const>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::DenseBase<Eigen::CwiseBinaryOp<Eigen::internal::scalar_cmp_op<int, (Eigen::internal::ComparisonName)0>, Eigen::Array<int, -1, 1, 0, -1, 1> const, Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<int>, Eigen::Array<int, -1, 1, 0, -1, 1> > const> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+#endif
+
 #endif

include/igl/flip_avoiding_line_search.cpp

@@ -72,9 +72,19 @@ namespace igl
           return INFINITY;
         }
 
-        double delta = sqrt(delta_in);
-        t1 = (-b + delta) / (2 * a);
-        t2 = (-b - delta) / (2 * a);
+        double delta = sqrt(delta_in); // delta >= 0
+        if(b >= 0) // avoid subtracting two similar numbers
+        {
+          double bd = - b - delta;
+          t1 = 2 * c / bd;
+          t2 = bd / (2 * a);
+        }
+        else
+        {
+          double bd = - b + delta;
+          t1 = bd / (2 * a);
+          t2 = (2 * c) / bd;
+        }
 
         assert (std::isfinite(t1));
         assert (std::isfinite(t2));
@@ -243,7 +253,7 @@ namespace igl
           return (res[0] >= 0) ? res[0]:INFINITY;
         case 2:
         {
-          double max_root = max(res[0],res[1]); double min_root = min(res[0],res[1]);
+          double max_root = std::max(res[0],res[1]); double min_root = std::min(res[0],res[1]);
           if (min_root > 0) return min_root;
           if (max_root > 0) return max_root;
           return INFINITY;
@@ -273,7 +283,7 @@ namespace igl
         for (int f = 0; f < F.rows(); f++)
         {
           double min_positive_root = get_min_pos_root_2D(uv,F,d,f);
-          max_step = min(max_step, min_positive_root);
+          max_step = std::min(max_step, min_positive_root);
         }
       }
       else
@@ -281,7 +291,7 @@ namespace igl
         for (int f = 0; f < F.rows(); f++)
         {
           double min_positive_root = get_min_pos_root_3D(uv,F,d,f);
-          max_step = min(max_step, min_positive_root);
+          max_step = std::min(max_step, min_positive_root);
         }
       }
       return max_step;
@@ -300,7 +310,7 @@ IGL_INLINE double igl::flip_avoiding_line_search(
   Eigen::MatrixXd d = dst_v - cur_v;
 
   double min_step_to_singularity = igl::flip_avoiding::compute_max_step_from_singularities(cur_v,F,d);
-  double max_step_size = min(1., min_step_to_singularity*0.8);
+  double max_step_size = std::min(1., min_step_to_singularity*0.8);
 
   return igl::line_search(cur_v,d,max_step_size, energy, cur_energy);
 }

include/igl/histc.cpp

@@ -100,10 +100,14 @@ IGL_INLINE void igl::histc(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
-template void igl::histc<Eigen::CwiseNullaryOp<Eigen::internal::linspaced_op<int, true>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::CwiseNullaryOp<Eigen::internal::linspaced_op<int, true>, Eigen::Matrix<int, -1, 1, 0, -1, 1> > > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::histc<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 template void igl::histc<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template void igl::histc<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::histc<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::histc<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+#if EIGEN_VERSION_AT_LEAST(3,3,0)
+#else
+template void igl::histc<Eigen::CwiseNullaryOp<Eigen::internal::linspaced_op<int, true>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::CwiseNullaryOp<Eigen::internal::linspaced_op<int, true>, Eigen::Matrix<int, -1, 1, 0, -1, 1> > > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+#endif
+
 #endif

include/igl/is_edge_manifold.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 "is_edge_manifold.h"
-#include "all_edges.h"
+#include "oriented_facets.h"
 #include "unique_simplices.h"
 
 #include <algorithm>
@@ -19,7 +19,7 @@ template <
   typename DerivedEMAP,
   typename DerivedBE>
 IGL_INLINE bool igl::is_edge_manifold(
-  const Eigen::PlainObjectBase<DerivedF>& F,
+  const Eigen::MatrixBase<DerivedF>& F,
   Eigen::PlainObjectBase<DerivedBF>& BF,
   Eigen::PlainObjectBase<DerivedE>& E,
   Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
@@ -30,7 +30,7 @@ IGL_INLINE bool igl::is_edge_manifold(
   typedef Matrix<typename DerivedF::Scalar,Dynamic,1> VectorXF;
   typedef Matrix<typename DerivedF::Scalar,Dynamic,2> MatrixXF2;
   MatrixXF2 allE;
-  all_edges(F,allE);
+  oriented_facets(F,allE);
   // Find unique undirected edges and mapping
   VectorXF _;
   unique_simplices(allE,E,_,EMAP);
@@ -54,7 +54,7 @@ IGL_INLINE bool igl::is_edge_manifold(
 
 template <typename DerivedF>
 IGL_INLINE bool igl::is_edge_manifold(
-  const Eigen::PlainObjectBase<DerivedF>& F)
+  const Eigen::MatrixBase<DerivedF>& F)
 {
   // TODO: It's bothersome that this is not calling/reusing the code from the
   // overload above. This could result in disagreement.
@@ -96,9 +96,9 @@ IGL_INLINE bool igl::is_edge_manifold(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
-template bool igl::is_edge_manifold<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&);
+template bool igl::is_edge_manifold<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&);
 // generated by autoexplicit.sh
 //template bool igl::is_edge_manifold<double>(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&);
-template bool igl::is_edge_manifold<Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&);
-template bool igl::is_edge_manifold<Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&);
+template bool igl::is_edge_manifold<Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&);
+template bool igl::is_edge_manifold<Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&);
 #endif

include/igl/is_edge_manifold.h

@@ -31,14 +31,14 @@ namespace igl
     typename DerivedEMAP,
     typename DerivedBE>
   IGL_INLINE bool is_edge_manifold(
-    const Eigen::PlainObjectBase<DerivedF>& F,
+    const Eigen::MatrixBase<DerivedF>& F,
     Eigen::PlainObjectBase<DerivedBF>& BF,
     Eigen::PlainObjectBase<DerivedE>& E,
     Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
     Eigen::PlainObjectBase<DerivedBE>& BE);
   template <typename DerivedF>
   IGL_INLINE bool is_edge_manifold(
-    const Eigen::PlainObjectBase<DerivedF>& F);
+    const Eigen::MatrixBase<DerivedF>& F);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/is_symmetric.cpp

@@ -34,13 +34,7 @@ IGL_INLINE bool igl::is_symmetric(
     return false;
   }
   assert(A.size() != 0);
-  const typename Eigen::PlainObjectBase<DerivedA>& AT = A.transpose();
-  const typename Eigen::PlainObjectBase<DerivedA>& AmAT = A-AT;
-  //// Eigen screws up something with LLT if you try to do
-  //SparseMatrix<T> AmAT = A-A.transpose();
-  //// Eigen crashes at runtime if you try to do
-  // return (A-A.transpose()).nonZeros() == 0;
-  return AmAT.nonZeros() == 0;
+  return (A-A.transpose()).eval().nonZeros() == 0;
 }
 
 template <typename AType, typename epsilonT>

include/igl/line_segment_in_rectangle.cpp

@@ -94,10 +94,10 @@ IGL_INLINE bool igl::line_segment_in_rectangle(
 
     return true;
   };
-  const double minX = min(A(0),B(0));
-  const double minY = min(A(1),B(1));
-  const double maxX = max(A(0),B(0));
-  const double maxY = max(A(1),B(1));
+  const double minX = std::min(A(0),B(0));
+  const double minY = std::min(A(1),B(1));
+  const double maxX = std::max(A(0),B(0));
+  const double maxY = std::max(A(1),B(1));
   bool ret = SegmentIntersectRectangle(minX,minY,maxX,maxY,s(0),s(1),d(0),d(1));
   return ret;
 }

include/igl/min_quad_with_fixed.cpp

@@ -406,8 +406,8 @@ IGL_INLINE bool igl::min_quad_with_fixed_solve(
   }
   // number of columns to solve
   int cols = Y.cols();
-  assert(B.cols() == 1);
-  assert(Beq.size() == 0 || Beq.cols() == 1);
+  assert(B.cols() == 1 || B.cols() == cols);
+  assert(Beq.size() == 0 || Beq.cols() == 1 || Beq.cols() == cols);
 
   // resize output
   Z.resize(data.n,cols);
@@ -425,25 +425,19 @@ IGL_INLINE bool igl::min_quad_with_fixed_solve(
     // number of lagrange multipliers aka linear equality constraints
     int neq = data.lagrange.size();
     // append lagrange multiplier rhs's
-    VectorXT BBeq(B.size() + Beq.size());
-    // Would like to do:
-    // BBeq << B, (Beq*-2.0);
-    // but Eigen can't handle empty vectors in comma initialization
-    // https://forum.kde.org/viewtopic.php?f=74&t=107974&p=364947#p364947
+    MatrixXT BBeq(B.rows() + Beq.rows(),cols);
     if(B.size() > 0)
     {
-      BBeq << B;
+      BBeq.topLeftCorner(B.rows(),cols) = B.replicate(1,B.cols()==cols?1:cols);
     }
     if(Beq.size() > 0)
     {
-      BBeq << Beq*-2.;
+      BBeq.bottomLeftCorner(Beq.rows(),cols) = -2.0*Beq.replicate(1,Beq.cols()==cols?1:cols);
     }
 
     // Build right hand side
-    VectorXT BBequl;
-    igl::slice(BBeq,data.unknown_lagrange,BBequl);
     MatrixXT BBequlcols;
-    igl::repmat(BBequl,1,cols,BBequlcols);
+    igl::slice(BBeq,data.unknown_lagrange,1,BBequlcols);
     MatrixXT NB;
     if(kr == 0)
     {
@@ -486,26 +480,26 @@ IGL_INLINE bool igl::min_quad_with_fixed_solve(
   }else
   {
     assert(data.solver_type == min_quad_with_fixed_data<T>::QR_LLT);
-    DerivedBeq eff_Beq;
+    MatrixXT eff_Beq;
     // Adjust Aeq rhs to include known parts
     eff_Beq =
       //data.AeqTQR.colsPermutation().transpose() * (-data.Aeqk * Y + Beq);
-      data.AeqTET * (-data.Aeqk * Y + Beq);
+      data.AeqTET * (-data.Aeqk * Y + Beq.replicate(1,Beq.cols()==cols?1:cols));
     // Where did this -0.5 come from? Probably the same place as above.
-    DerivedB Bu;
-    slice(B,data.unknown,Bu);
-    DerivedB NB;
-    NB = -0.5*(Bu + data.preY * Y);
+    MatrixXT Bu;
+    slice(B,data.unknown,1,Bu);
+    MatrixXT NB;
+    NB = -0.5*(Bu.replicate(1,B.cols()==cols?1:cols) + data.preY * Y);
     // Trim eff_Beq
     const int nc = data.AeqTQR.rank();
     const int neq = Beq.rows();
-    eff_Beq = eff_Beq.topLeftCorner(nc,1).eval();
+    eff_Beq = eff_Beq.topLeftCorner(nc,cols).eval();
     data.AeqTR1T.template triangularView<Lower>().solveInPlace(eff_Beq);
     // Now eff_Beq = (data.AeqTR1T \ (data.AeqTET * (-data.Aeqk * Y + Beq)))
-    DerivedB lambda_0;
+    MatrixXT lambda_0;
     lambda_0 = data.AeqTQ1 * eff_Beq;
     //cout<<matlab_format(lambda_0,"lambda_0")<<endl;
-    DerivedB QRB;
+    MatrixXT QRB;
     QRB = -data.AeqTQ2T * (data.Auu * lambda_0) + data.AeqTQ2T * NB;
     Derivedsol lambda;
     lambda = data.llt.solve(QRB);
@@ -519,8 +513,8 @@ IGL_INLINE bool igl::min_quad_with_fixed_solve(
       temp1 = (data.AeqTQ1T * NB - data.AeqTQ1T * data.Auu * solu);
       data.AeqTR1.template triangularView<Upper>().solveInPlace(temp1);
       //cout<<matlab_format(temp1,"temp1")<<endl;
-      temp2 = Derivedsol::Zero(neq,1);
-      temp2.topLeftCorner(nc,1) = temp1;
+      temp2 = Derivedsol::Zero(neq,cols);
+      temp2.topLeftCorner(nc,cols) = temp1;
       //solLambda = data.AeqTQR.colsPermutation() * temp2;
       solLambda = data.AeqTE * temp2;
     }

include/igl/min_quad_with_fixed.h

@@ -26,12 +26,14 @@ namespace igl
   // they're not then resulting probably will no longer be sparse: it will be
   // slow.
   //
-  // MIN_QUAD_WITH_FIXED Minimize quadratic energy 
+  // MIN_QUAD_WITH_FIXED Minimize a quadratic energy of the form
   //
-  // 0.5*Z'*A*Z + Z'*B + C with
+  // trace( 0.5*Z'*A*Z + Z'*B + constant )
   //
-  // constraints that Z(known) = Y, optionally also subject to the constraints
-  // Aeq*Z = Beq
+  // subject to
+  //
+  //   Z(known,:) = Y, and
+  //   Aeq*Z = Beq
   //
   // Templates:
   //   T  should be a eigen matrix primitive type like int or double
@@ -57,7 +59,6 @@ namespace igl
     const bool pd,
     min_quad_with_fixed_data<T> & data
     );
-
   // Solves a system previously factored using min_quad_with_fixed_precompute
   //
   // Template:
@@ -66,12 +67,12 @@ namespace igl
   //   DerivedZ  type of Z (e.g. derived from VectorXd or MatrixXd)
   // Inputs:
   //   data  factorization struct with all necessary precomputation to solve
-  //   B  n by 1 column of linear coefficients
-  //   Y  b by 1 list of constant fixed values
-  //   Beq  m by 1 list of linear equality constraint constant values
+  //   B  n by k column of linear coefficients
+  //   Y  b by k list of constant fixed values
+  //   Beq  m by k list of linear equality constraint constant values
   // Outputs:
-  //   Z  n by cols solution
-  //   sol  #unknowns+#lagrange by cols solution to linear system
+  //   Z  n by k solution
+  //   sol  #unknowns+#lagrange by k solution to linear system
   // Returns true on success, false on error
   template <
     typename T,

include/igl/mosek/mosek_quadprog.cpp

@@ -282,6 +282,8 @@ IGL_INLINE bool igl::mosek::mosek_quadprog(
   using namespace Eigen;
   using namespace std;
 
+  typedef int Index;
+  typedef double Scalar;
   // Q should be square
   assert(Q.rows() == Q.cols());
   // Q should be symmetric
@@ -289,36 +291,36 @@ IGL_INLINE bool igl::mosek::mosek_quadprog(
   assert( (Q-Q.transpose()).sum() < FLOAT_EPS);
 #endif
   // Only keep lower triangular part of Q
-  SparseMatrix<double> QL;
+  SparseMatrix<Scalar> QL;
   //QL = Q.template triangularView<Lower>();
   QL = Q.triangularView<Lower>();
   VectorXi Qi,Qj;
   VectorXd Qv;
   find(QL,Qi,Qj,Qv);
-  vector<int> vQi = matrix_to_list(Qi);
-  vector<int> vQj = matrix_to_list(Qj);
-  vector<double> vQv = matrix_to_list(Qv);
+  vector<Index> vQi = matrix_to_list(Qi);
+  vector<Index> vQj = matrix_to_list(Qj);
+  vector<Scalar> vQv = matrix_to_list(Qv);
 
   // Convert linear term
-  vector<double> vc = matrix_to_list(c);
+  vector<Scalar> vc = matrix_to_list(c);
 
   assert(lc.size() == A.rows());
   assert(uc.size() == A.rows());
   // Convert A to harwell boeing format
-  vector<double> vAv;
-  vector<int> vAr,vAc;
-  int nr;
+  vector<Scalar> vAv;
+  vector<Index> vAr,vAc;
+  Index nr;
   harwell_boeing(A,nr,vAv,vAr,vAc);
 
   assert(lx.size() == Q.rows());
   assert(ux.size() == Q.rows());
-  vector<double> vlc = matrix_to_list(lc);
-  vector<double> vuc = matrix_to_list(uc);
-  vector<double> vlx = matrix_to_list(lx);
-  vector<double> vux = matrix_to_list(ux);
+  vector<Scalar> vlc = matrix_to_list(lc);
+  vector<Scalar> vuc = matrix_to_list(uc);
+  vector<Scalar> vlx = matrix_to_list(lx);
+  vector<Scalar> vux = matrix_to_list(ux);
 
-  vector<double> vx;
-  bool ret = mosek_quadprog(
+  vector<Scalar> vx;
+  bool ret = mosek_quadprog<Index,Scalar>(
     Q.rows(),vQi,vQj,vQv,
     vc,
     cf,

include/igl/nchoosek.cpp

@@ -68,6 +68,10 @@ IGL_INLINE void igl::nchoosek(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
-template void igl::nchoosek<Eigen::CwiseNullaryOp<Eigen::internal::linspaced_op<int, true>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::CwiseNullaryOp<Eigen::internal::linspaced_op<int, true>, Eigen::Matrix<int, -1, 1, 0, -1, 1> > > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template void igl::nchoosek<Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+#if EIGEN_VERSION_AT_LEAST(3,3,0)
+#else
+template void igl::nchoosek<Eigen::CwiseNullaryOp<Eigen::internal::linspaced_op<int, true>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::CwiseNullaryOp<Eigen::internal::linspaced_op<int, true>, Eigen::Matrix<int, -1, 1, 0, -1, 1> > > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+#endif
+
 #endif

include/igl/next_filename.h

@@ -12,12 +12,12 @@
 namespace igl
 {
   // Find the file with the first filename of the form
-  // "prefix-%0[zeros]dsuffix"
+  // "prefix%0[zeros]dsuffix"
   // 
   // Inputs:
   //   prefix  path to containing dir and filename prefix
   //   zeros number of leading zeros as if digit printed with printf
-  //   suffix  suffix of filename and extension
+  //   suffix  suffix of filename and extension (should included dot)
   // Outputs:
   //   next  path to next file
   // Returns true if found, false if exceeding range in zeros

include/igl/oriented_facets.cpp

@@ -0,0 +1,55 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2017 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 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "oriented_facets.h"
+
+template <typename DerivedF, typename DerivedE>
+IGL_INLINE void igl::oriented_facets(
+  const Eigen::MatrixBase<DerivedF> & F,
+  Eigen::PlainObjectBase<DerivedE> & E)
+{
+  E.resize(F.rows()*F.cols(),F.cols()-1);
+  typedef typename DerivedE::Scalar EScalar;
+  switch(F.cols())
+  {
+    case 4:
+      E.block(0*F.rows(),0,F.rows(),1) = F.col(1).template cast<EScalar>();
+      E.block(0*F.rows(),1,F.rows(),1) = F.col(3).template cast<EScalar>();
+      E.block(0*F.rows(),2,F.rows(),1) = F.col(2).template cast<EScalar>();
+
+      E.block(1*F.rows(),0,F.rows(),1) = F.col(0).template cast<EScalar>();
+      E.block(1*F.rows(),1,F.rows(),1) = F.col(2).template cast<EScalar>();
+      E.block(1*F.rows(),2,F.rows(),1) = F.col(3).template cast<EScalar>();
+
+      E.block(2*F.rows(),0,F.rows(),1) = F.col(0).template cast<EScalar>();
+      E.block(2*F.rows(),1,F.rows(),1) = F.col(3).template cast<EScalar>();
+      E.block(2*F.rows(),2,F.rows(),1) = F.col(1).template cast<EScalar>();
+
+      E.block(3*F.rows(),0,F.rows(),1) = F.col(0).template cast<EScalar>();
+      E.block(3*F.rows(),1,F.rows(),1) = F.col(1).template cast<EScalar>();
+      E.block(3*F.rows(),2,F.rows(),1) = F.col(2).template cast<EScalar>();
+      return;
+    case 3:
+      E.block(0*F.rows(),0,F.rows(),1) = F.col(1).template cast<EScalar>();
+      E.block(0*F.rows(),1,F.rows(),1) = F.col(2).template cast<EScalar>();
+      E.block(1*F.rows(),0,F.rows(),1) = F.col(2).template cast<EScalar>();
+      E.block(1*F.rows(),1,F.rows(),1) = F.col(0).template cast<EScalar>();
+      E.block(2*F.rows(),0,F.rows(),1) = F.col(0).template cast<EScalar>();
+      E.block(2*F.rows(),1,F.rows(),1) = F.col(1).template cast<EScalar>();
+      return;
+  }
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template void igl::oriented_facets<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::oriented_facets<Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 2, 0, -1, 2> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&);
+template void igl::oriented_facets<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 2, 0, -1, 2> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&);
+template void igl::oriented_facets<Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::oriented_facets<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 2, 0, -1, 2> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> >&);
+#endif
+

include/igl/oriented_facets.h

@@ -0,0 +1,41 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2017 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 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_ORIENTED_FACETS_H
+#define IGL_ORIENTED_FACETS_H
+#include "igl_inline.h"
+#include <Eigen/Dense>
+namespace igl
+{
+  // ORIENTED_FACETS Determines all "directed
+  // [facets](https://en.wikipedia.org/wiki/Simplex#Elements)" of a given set
+  // of simplicial elements. For a manifold triangle mesh, this computes all
+  // half-edges. For a manifold tetrahedral mesh, this computes all half-faces.
+  //
+  // Inputs:
+  //   F  #F by simplex_size list of simplices
+  // Outputs:
+  //   E  #E by simplex_size-1  list of facets
+  //
+  // Note: this is not the same as igl::edges because this includes every
+  // directed edge including repeats (meaning interior edges on a surface will
+  // show up once for each direction and non-manifold edges may appear more than
+  // once for each direction).
+  //
+  // Note: This replaces the deprecated `all_edges` function
+  template <typename DerivedF, typename DerivedE>
+  IGL_INLINE void oriented_facets(
+    const Eigen::MatrixBase<DerivedF> & F,
+    Eigen::PlainObjectBase<DerivedE> & E);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "oriented_facets.cpp"
+#endif
+
+#endif
+

include/igl/per_edge_normals.cpp

@@ -5,7 +5,7 @@
 // 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 "all_edges.h"
+#include "oriented_facets.h"
 #include "doublearea.h"
 #include "per_edge_normals.h"
 #include "get_seconds.h"
@@ -37,7 +37,7 @@ IGL_INLINE void igl::per_edge_normals(
   const int m = F.rows();
   // All occurances of directed edges
   MatrixXi allE;
-  all_edges(F,allE);
+  oriented_facets(F,allE);
   // Find unique undirected edges and mapping
   VectorXi _;
   unique_simplices(allE,E,_,EMAP);

include/igl/pinv.cpp

@@ -0,0 +1,35 @@
+#include "pinv.h"
+#include <limits>
+#include <cmath>
+
+template <typename DerivedA, typename DerivedX>
+void igl::pinv(
+  const Eigen::MatrixBase<DerivedA> & A,
+  typename DerivedA::Scalar tol,
+  Eigen::PlainObjectBase<DerivedX> & X)
+{
+  Eigen::JacobiSVD<DerivedA> svd(A, Eigen::ComputeFullU | Eigen::ComputeFullV );
+  typedef typename DerivedA::Scalar Scalar;
+  const Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic> & U = svd.matrixU();
+  const Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic> & V = svd.matrixV();
+  const Eigen::Matrix<Scalar,Eigen::Dynamic,1> & S = svd.singularValues();
+  if(tol < 0)
+  {
+    const Scalar smax = S.array().abs().maxCoeff();
+    tol = 
+      (Scalar)(std::max(A.rows(),A.cols())) *
+      (smax-std::nextafter(smax,std::numeric_limits<Scalar>::epsilon()));
+  }
+  const int rank = (S.array()>0).count();
+  X = (V.leftCols(rank).array().rowwise() * 
+      (1.0/S.head(rank).array()).transpose()).matrix()*
+    U.leftCols(rank).transpose();
+}
+
+template <typename DerivedA, typename DerivedX>
+void igl::pinv(
+  const Eigen::MatrixBase<DerivedA> & A,
+  Eigen::PlainObjectBase<DerivedX> & X)
+{
+  return pinv(A,-1,X);
+}

include/igl/pinv.h

@@ -0,0 +1,29 @@
+#ifndef IGL_PINV_H
+#define IGL_PINV_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // Compute the Moore-Penrose pseudoinverse
+  //
+  // Inputs:
+  //   A  m by n matrix 
+  //   tol  tolerance (if negative then default is used)
+  // Outputs:
+  //   X  n by m matrix so that A*X*A = A and X*A*X = X and A*X = (A*X)' and
+  //     (X*A) = (X*A)'
+  template <typename DerivedA, typename DerivedX>
+  void pinv(
+    const Eigen::MatrixBase<DerivedA> & A,
+    typename DerivedA::Scalar tol,
+    Eigen::PlainObjectBase<DerivedX> & X);
+  // Wrapper using default tol
+  template <typename DerivedA, typename DerivedX>
+  void pinv(
+    const Eigen::MatrixBase<DerivedA> & A,
+    Eigen::PlainObjectBase<DerivedX> & X);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "pinv.cpp"
+#endif
+#endif

include/igl/readOFF.cpp

@@ -67,7 +67,7 @@ IGL_INLINE bool igl::readOFF(
   while(still_comments)
   {
     fgets(line,1000,off_file);
-    still_comments = line[0] == '#';
+    still_comments = (line[0] == '#' || line[0] == '\n');
   }
   sscanf(line,"%d %d %d",&number_of_vertices,&number_of_faces,&number_of_edges);
   V.resize(number_of_vertices);

include/igl/remove_unreferenced.cpp

@@ -45,7 +45,7 @@ IGL_INLINE void igl::remove_unreferenced(
   const size_t n = V.rows();
   remove_unreferenced(n,F,I,J);
   NF = F;
-  for_each(NF.data(),NF.data()+NF.size(),
+  std::for_each(NF.data(),NF.data()+NF.size(),
     [&I](typename DerivedNF::Scalar & a){a=I(a);});
   slice(V,J,1,NV);
 }

include/igl/setdiff.cpp

@@ -42,7 +42,8 @@ IGL_INLINE void igl::setdiff(
       }
     }
     assert(k == C.size());
-    IA = igl::LinSpaced<Eigen::Matrix<typename DerivedIA::Scalar,Eigen::Dynamic,1> >(
+    // Have to use << instead of = becasue Eigen's PlainObjectBase is annoying
+    IA << igl::LinSpaced<Eigen::Matrix<typename DerivedIA::Scalar,Eigen::Dynamic,1> >(
       C.size(),0,C.size()-1);
   }
 

include/igl/slice.cpp

@@ -248,6 +248,12 @@ IGL_INLINE DerivedX igl::slice(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template void igl::slice<Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int, Eigen::Matrix<double, -1, -1, 0, -1, -1>&);
+// generated by autoexplicit.sh
+template void igl::slice<Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int, Eigen::Matrix<double, -1, -1, 0, -1, -1>&);
+// generated by autoexplicit.sh
+template void igl::slice<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > >(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
 template void igl::slice<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 1, -1, -1> >(Eigen::Matrix<double, -1, -1, 1, -1, -1> const&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int, Eigen::Matrix<double, -1, -1, 1, -1, -1>&);
 // generated by autoexplicit.sh
 template void igl::slice<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 1, -1, -1> >(Eigen::DenseBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> >&);

include/igl/sort_angles.cpp

@@ -18,7 +18,8 @@ IGL_INLINE void igl::sort_angles(
     assert(num_cols >= 2);
 
     R.resize(num_rows);
-    R = igl::LinSpaced<
+    // Have to use << instead of = becasue Eigen's PlainObjectBase is annoying
+    R << igl::LinSpaced<
       Eigen::Matrix<typename DerivedR::Scalar, Eigen::Dynamic, 1> >
       (num_rows, 0, num_rows-1);
 

include/igl/sparse.cpp

@@ -110,10 +110,13 @@ IGL_INLINE Eigen::SparseMatrix<typename DerivedD::Scalar > igl::sparse(
 // generated by autoexplicit.sh
 template void igl::sparse<Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<bool>, Eigen::Array<bool, -1, 2, 0, -1, 2> >, bool>(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<bool>, Eigen::Array<bool, -1, 2, 0, -1, 2> > const&, unsigned long, unsigned long, Eigen::SparseMatrix<bool, 0, int>&);
 // generated by autoexplicit.sh
-template void igl::sparse<Eigen::CwiseNullaryOp<Eigen::internal::linspaced_op<int, true>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<bool>, Eigen::Array<bool, -1, 1, 0, -1, 1> >, bool>(Eigen::CwiseNullaryOp<Eigen::internal::linspaced_op<int, true>, Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<bool>, Eigen::Array<bool, -1, 1, 0, -1, 1> > const&, unsigned long, unsigned long, Eigen::SparseMatrix<bool, 0, int>&);
 // generated by autoexplicit.sh
 template void igl::sparse<Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >, Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<bool>, Eigen::Array<bool, -1, 3, 0, -1, 3> >, bool>(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<bool>, Eigen::Array<bool, -1, 3, 0, -1, 3> > const&, unsigned long, unsigned long, Eigen::SparseMatrix<bool, 0, int>&);
 template void igl::sparse<Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, std::complex<double> >(Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::Matrix<double, -1, 1, 0, -1, 1> const&, size_t, size_t, Eigen::SparseMatrix<std::complex<double>, 0, int>&);
 template void igl::sparse<Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, double>(Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::Matrix<double, -1, 1, 0, -1, 1> const&, Eigen::SparseMatrix<double, 0, int>&);
 template void igl::sparse<Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, double>(Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::Matrix<double, -1, 1, 0, -1, 1> const&, size_t, size_t, Eigen::SparseMatrix<double, 0, int>&);
+#if EIGEN_VERSION_AT_LEAST(3,3,0)
+#else
+template void igl::sparse<Eigen::CwiseNullaryOp<Eigen::internal::linspaced_op<int, true>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<bool>, Eigen::Array<bool, -1, 1, 0, -1, 1> >, bool>(Eigen::CwiseNullaryOp<Eigen::internal::linspaced_op<int, true>, Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<bool>, Eigen::Array<bool, -1, 1, 0, -1, 1> > const&, unsigned long, unsigned long, Eigen::SparseMatrix<bool, 0, int>&);
+#endif
 #endif

include/igl/sum.cpp

@@ -57,6 +57,8 @@ IGL_INLINE void igl::sum(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template void igl::sum<double, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::SparseMatrix<double, 0, int> const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+// generated by autoexplicit.sh
 template void igl::sum<bool, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::SparseMatrix<bool, 0, int> const&, int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::sum<double>(Eigen::SparseMatrix<double, 0, int> const&, int, Eigen::SparseVector<double, 0, int>&);
 #endif

include/igl/swept_volume_signed_distance.cpp

@@ -16,6 +16,7 @@
 #include "per_edge_normals.h"
 #include <Eigen/Geometry>
 #include <cmath>
+#include <algorithm>
 
 IGL_INLINE void igl::swept_volume_signed_distance(
   const Eigen::MatrixXd & V,
@@ -80,7 +81,7 @@ IGL_INLINE void igl::swept_volume_signed_distance(
         pseudonormal_test(V,F,FN,VN,EN,EMAP,gv,i,c,s,n);
         if(S(g) == S(g))
         {
-          S(g) = min(S(g),s*sqrt(sqrd));
+          S(g) = std::min(S(g),s*sqrt(sqrd));
         }else
         {
           S(g) = s*sqrt(sqrd);
@@ -96,7 +97,7 @@ IGL_INLINE void igl::swept_volume_signed_distance(
   {
 #ifndef NDEBUG
     // Check for nans
-    for_each(S.data(),S.data()+S.size(),[](const double s){assert(s==s);});
+    std::for_each(S.data(),S.data()+S.size(),[](const double s){assert(s==s);});
 #endif
   }
 }

include/igl/triangle_triangle_adjacency.cpp

@@ -6,8 +6,6 @@
 // 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 "triangle_triangle_adjacency.h"
-#include "all_edges.h"
-#include "unique_simplices.h"
 #include "parallel_for.h"
 #include "unique_edge_map.h"
 #include <algorithm>

include/igl/triangle_triangle_adjacency.h

@@ -90,11 +90,12 @@ namespace igl
       std::vector<std::vector<std::vector<TTIndex> > > & TT,
       std::vector<std::vector<std::vector<TTiIndex> > > & TTi);
   // Inputs:
-  //   E  #F*3 by 2 list of all of directed edges in order (see `all_edges`)
+  //   E  #F*3 by 2 list of all of directed edges in order (see
+  //     `oriented_facets`)
   //   EMAP #F*3 list of indices into uE, mapping each directed edge to unique
   //     undirected edge
   //   uE2E  #uE list of lists of indices into E of coexisting edges
-  // See also: unique_edge_map, all_edges
+  // See also: unique_edge_map, oriented_facets
   template <
     typename DerivedE, 
     typename DerivedEMAP,

include/igl/unique_edge_map.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 "unique_edge_map.h"
-#include "all_edges.h"
+#include "oriented_facets.h"
 #include "unique_simplices.h"
 #include <cassert>
 #include <algorithm>
@@ -26,7 +26,7 @@ IGL_INLINE void igl::unique_edge_map(
   using namespace Eigen;
   using namespace std;
   // All occurances of directed edges
-  all_edges(F,E);
+  oriented_facets(F,E);
   const size_t ne = E.rows();
   // This is 2x faster to create than a map from pairs to lists of edges and 5x
   // faster to access (actually access is probably assympotically faster O(1)

include/igl/unproject.cpp

@@ -17,10 +17,10 @@ template <
   typename Derivedviewport,
   typename Derivedscene>
 IGL_INLINE void igl::unproject(
-  const Eigen::PlainObjectBase<Derivedwin>&  win,
-  const Eigen::PlainObjectBase<Derivedmodel>& model,
-  const Eigen::PlainObjectBase<Derivedproj>& proj,
-  const Eigen::PlainObjectBase<Derivedviewport>&  viewport,
+  const Eigen::MatrixBase<Derivedwin>&  win,
+  const Eigen::MatrixBase<Derivedmodel>& model,
+  const Eigen::MatrixBase<Derivedproj>& proj,
+  const Eigen::MatrixBase<Derivedviewport>&  viewport,
   Eigen::PlainObjectBase<Derivedscene> & scene)
 {
   if(win.cols() != 3)
@@ -69,6 +69,6 @@ IGL_INLINE Eigen::Matrix<Scalar,3,1> igl::unproject(
 #ifdef IGL_STATIC_LIBRARY
 template Eigen::Matrix<float, 3, 1, 0, 3, 1> igl::unproject<float>(Eigen::Matrix<float, 3, 1, 0, 3, 1> const&, Eigen::Matrix<float, 4, 4, 0, 4, 4> const&, Eigen::Matrix<float, 4, 4, 0, 4, 4> const&, Eigen::Matrix<float, 4, 1, 0, 4, 1> const&);
 template Eigen::Matrix<double, 3, 1, 0, 3, 1> igl::unproject<double>(Eigen::Matrix<double, 3, 1, 0, 3, 1> const&, Eigen::Matrix<double, 4, 4, 0, 4, 4> const&, Eigen::Matrix<double, 4, 4, 0, 4, 4> const&, Eigen::Matrix<double, 4, 1, 0, 4, 1> const&);
-template void igl::unproject<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<float, 4, 4, 0, 4, 4>, Eigen::Matrix<float, 4, 4, 0, 4, 4>, Eigen::Matrix<float, 4, 1, 0, 4, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, 4, 4, 0, 4, 4> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, 4, 4, 0, 4, 4> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, 4, 1, 0, 4, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
-template void igl::unproject<Eigen::Matrix<float, 3, 1, 0, 3, 1>, Eigen::Matrix<float, 4, 4, 0, 4, 4>, Eigen::Matrix<float, 4, 4, 0, 4, 4>, Eigen::Matrix<float, 4, 1, 0, 4, 1>, Eigen::Matrix<float, 3, 1, 0, 3, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<float, 3, 1, 0, 3, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, 4, 4, 0, 4, 4> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, 4, 4, 0, 4, 4> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, 4, 1, 0, 4, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, 3, 1, 0, 3, 1> >&);
+template void igl::unproject<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<float, 4, 4, 0, 4, 4>, Eigen::Matrix<float, 4, 4, 0, 4, 4>, Eigen::Matrix<float, 4, 1, 0, 4, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<float, 4, 4, 0, 4, 4> > const&, Eigen::MatrixBase<Eigen::Matrix<float, 4, 4, 0, 4, 4> > const&, Eigen::MatrixBase<Eigen::Matrix<float, 4, 1, 0, 4, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::unproject<Eigen::Matrix<float, 3, 1, 0, 3, 1>, Eigen::Matrix<float, 4, 4, 0, 4, 4>, Eigen::Matrix<float, 4, 4, 0, 4, 4>, Eigen::Matrix<float, 4, 1, 0, 4, 1>, Eigen::Matrix<float, 3, 1, 0, 3, 1> >(Eigen::MatrixBase<Eigen::Matrix<float, 3, 1, 0, 3, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<float, 4, 4, 0, 4, 4> > const&, Eigen::MatrixBase<Eigen::Matrix<float, 4, 4, 0, 4, 4> > const&, Eigen::MatrixBase<Eigen::Matrix<float, 4, 1, 0, 4, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, 3, 1, 0, 3, 1> >&);
 #endif

include/igl/unproject.h

@@ -27,10 +27,10 @@ namespace igl
     typename Derivedviewport,
     typename Derivedscene>
   IGL_INLINE void unproject(
-    const Eigen::PlainObjectBase<Derivedwin>&  win,
-    const Eigen::PlainObjectBase<Derivedmodel>& model,
-    const Eigen::PlainObjectBase<Derivedproj>& proj,
-    const Eigen::PlainObjectBase<Derivedviewport>&  viewport,
+    const Eigen::MatrixBase<Derivedwin>&  win,
+    const Eigen::MatrixBase<Derivedmodel>& model,
+    const Eigen::MatrixBase<Derivedproj>& proj,
+    const Eigen::MatrixBase<Derivedviewport>&  viewport,
     Eigen::PlainObjectBase<Derivedscene> & scene);
   template <typename Scalar>
   IGL_INLINE Eigen::Matrix<Scalar,3,1> unproject(

include/igl/vector_area_matrix.h

@@ -5,8 +5,8 @@
 // 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/.
-#ifndef IGL_AREAMATRIX_H
-#define IGL_AREAMATRIX_H
+#ifndef IGL_VECTOR_AREA_MATRIX_H
+#define IGL_VECTOR_AREA_MATRIX_H
 #include "igl_inline.h"
 
 #include <Eigen/Dense>

include/igl/volume.cpp

@@ -22,10 +22,11 @@ IGL_INLINE void igl::volume(
   vol.resize(m,1);
   for(int t = 0;t<m;t++)
   {
-    const RowVector3d & a = V.row(T(t,0));
-    const RowVector3d & b = V.row(T(t,1));
-    const RowVector3d & c = V.row(T(t,2));
-    const RowVector3d & d = V.row(T(t,3));
+    typedef Eigen::Matrix<typename DerivedV::Scalar,1,3> RowVector3S;
+    const RowVector3S & a = V.row(T(t,0));
+    const RowVector3S & b = V.row(T(t,1));
+    const RowVector3S & c = V.row(T(t,2));
+    const RowVector3S & d = V.row(T(t,3));
     vol(t) = -(a-d).dot((b-d).cross(c-d))/6.;
   }
 }
@@ -107,6 +108,8 @@ IGL_INLINE void igl::volume(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template void igl::volume<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+// generated by autoexplicit.sh
 template void igl::volume<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 // generated by autoexplicit.sh
 template void igl::volume<Eigen::Matrix<double, -1, 6, 0, -1, 6>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 6, 0, -1, 6> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);

shared/cmake/CMakeLists.txt

@@ -65,7 +65,9 @@ set(LIBIGL_EXTERNAL "${LIBIGL_ROOT}/external")
 set(NANOGUI_DIR "${LIBIGL_EXTERNAL}/nanogui")
 
 ### Eigen ###
-set(EIGEN_INCLUDE_DIR "${NANOGUI_DIR}/ext/eigen")
+if(NOT EIGEN_INCLUDE_DIR)
+  set(EIGEN_INCLUDE_DIR "${NANOGUI_DIR}/ext/eigen")
+endif()
 list(APPEND LIBIGL_INCLUDE_DIRS "${EIGEN_INCLUDE_DIR}")
 
 macro(CompileIGL_Module module_dir prefix)
@@ -152,7 +154,9 @@ endif()
 
 ### Compile the cgal parts ###
 if(LIBIGL_WITH_CGAL) # to be cleaned
-  find_package(CGAL REQUIRED)
+  # Core is needed for
+  # `Exact_predicates_exact_constructions_kernel_with_sqrt`
+  find_package(CGAL REQUIRED COMPONENTS Core)
   # set(Boost_USE_MULTITHREADED      ON)
   # set(Boost_USE_STATIC_LIBS      ON)
   #
@@ -334,6 +338,7 @@ if(LIBIGL_WITH_OPENGL)
     list(APPEND LIBIGL_INCLUDE_DIRS "${NANOGUI_DIR}/ext/glew/include")
     list(APPEND LIBIGL_OPENGL_EXTRA_LIBRARIES "glew")
   endif()
+  list(APPEND LIBIGL_EXTRA_LIBRARIES ${LIBIGL_OPENGL_EXTRA_LIBRARIES})
 
   list(APPEND LIBIGL_EXTRA_LIBRARIES ${LIBIGL_OPENGL_EXTRA_LIBRARIES})
 

shared/cmake/FindMATLAB.cmake

@@ -137,7 +137,7 @@ ELSE(WIN32)
     IF((NOT DEFINED MATLAB_ROOT) OR ("${MATLAB_ROOT}" STREQUAL ""))
 
     # Search for a version of Matlab available, starting from the most modern one to older versions
-      FOREACH(MATVER "R2016b" "R2016a" "R2015b" "R2015a" "R2014b" "R2014a" "R2014a" "R2013b" "R2013a" "R2012b" "R2012a" "R2011b" "R2011a" "R2010b" "R2010a" "R2009b" "R2009a" "R2008b")
+      FOREACH(MATVER "R2017b" "R2017a" "R2016b" "R2016a" "R2015b" "R2015a" "R2014b" "R2014a" "R2014a" "R2013b" "R2013a" "R2012b" "R2012a" "R2011b" "R2011a" "R2010b" "R2010a" "R2009b" "R2009a" "R2008b")
         IF((NOT DEFINED MATLAB_ROOT) OR ("${MATLAB_ROOT}" STREQUAL ""))
           IF(EXISTS /Applications/MATLAB_${MATVER}.app)
             SET(MATLAB_ROOT /Applications/MATLAB_${MATVER}.app)

tutorial/CMakeLists.txt

@@ -9,7 +9,7 @@ option(LIBIGL_WITH_VIEWER      "Use OpenGL viewer"  ON)
 option(LIBIGL_WITH_NANOGUI     "Use Nanogui menu"   OFF)
 
 ### libIGL options: choose your dependencies (by default everything is OFF, in this example we need the viewer) ###
-find_package(CGAL QUIET)
+find_package(CGAL QUIET COMPONENTS Core)
 option(LIBIGL_WITH_CGAL             "Use CGAL"           "${CGAL_FOUND}")
 option(LIBIGL_WITH_COMISO           "Use CoMiso"         ON)
 option(LIBIGL_WITH_CORK             "Use CORK"           OFF)
@@ -28,7 +28,7 @@ option(LIBIGL_WITH_XML              "Use XML"            ON)
 
 ### libIGL options: decide if you want to use the functionalities that depends on cgal
 if(LIBIGL_WITH_CGAL) # Do not remove or move this block, cgal strange build system fails without it
-  find_package(CGAL REQUIRED)
+  find_package(CGAL REQUIRED COMPONENTS Core)
   set(CGAL_DONT_OVERRIDE_CMAKE_FLAGS TRUE CACHE BOOL "CGAL's CMAKE Setup is super annoying ")
   include(${CGAL_USE_FILE})
 endif()

tutorial/tutorial.md.REMOVED.git-id

@@ -1 +1 @@
-7e9a0845653ff898c4654b6c35e9d326955283f8
+70fb600fc4ac9cf8080fc3dd5b25409140700b8a