Merge branch 'master' of https://github.com/libigl/libigl

Former-commit-id: eacb7c408a2d63f9bfdfb7835c88bf6ef97f781b

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/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/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/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/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/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/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/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> >&);