Unittests!! (#928)

* Migrate unit tests back to the main repo.

* Minor bug fixes and getting unit test migrated.


Former-commit-id: 823ef3e404ca32ec5e70c2d44bf0664b62631032

cmake/LibiglDownloadExternal.cmake

@@ -124,7 +124,7 @@ endfunction()
 function(igl_download_googletest)
 	igl_download_project(googletest
 		GIT_REPOSITORY https://github.com/google/googletest
-		GIT_TAG        v1.8.1
+		GIT_TAG        release-1.8.1
 	)
 endfunction()
 

include/igl/cut_to_disk.cpp

@@ -9,7 +9,7 @@ namespace igl {
   template <typename DerivedF, typename Index>
   void cut_to_disk(
     const Eigen::MatrixBase<DerivedF> &F,
-    std::vector<std::vector<Index> > &cuts) 
+    std::vector<std::vector<Index> > &cuts)
   {
     cuts.clear();
 
@@ -20,7 +20,7 @@ namespace igl {
 
     std::map<std::pair<Index, Index>, std::vector<Index> > edges;
     // build edges
-    
+
     for (Index i = 0; i < nfaces; i++)
     {
         for (int j = 0; j < 3; j++)
@@ -70,13 +70,13 @@ namespace igl {
             faceEdges(i, j) = edgeidx[e];
         }
     }
-    
+
     bool *deleted = new bool[nfaces];
     for (Index i = 0; i < nfaces; i++)
         deleted[i] = false;
 
     std::set<Index> deletededges;
-    
+
     // loop over faces
     for (Index face = 0; face < nfaces; face++)
     {
@@ -146,7 +146,7 @@ namespace igl {
         spinevertedges[edgeVerts(i, 0)].push_back(i);
         spinevertedges[edgeVerts(i, 1)].push_back(i);
     }
-    
+
     std::deque<Index> vertsProcess;
     std::map<Index, int> spinevertnbs;
     for (auto it : spinevertedges)
@@ -190,7 +190,7 @@ namespace igl {
         loopvertedges[edgeVerts(e, 0)].push_back(e);
         loopvertedges[edgeVerts(e, 1)].push_back(e);
     }
-    
+
     std::set<Index> usededges;
     for (Index e : loopedges)
     {
@@ -250,9 +250,9 @@ namespace igl {
                     {
                         cycleidx[nextvert] = cycleverts.size();
                         cycleverts.push_back(nextvert);
-                        cycleedges.push_back(nexte);                        
+                        cycleedges.push_back(nexte);
                     }
-                }                
+                }
                 curvert = nextvert;
                 cure = nexte;
             }
@@ -303,9 +303,9 @@ namespace igl {
                         {
                             cycleidx[nextvert] = cycleverts.size();
                             cycleverts.push_back(nextvert);
-                            cycleedges.push_back(nexte);                        
+                            cycleedges.push_back(nexte);
                         }
-                    }                
+                    }
                     curvert = nextvert;
                     cure = nexte;
                 }
@@ -321,10 +321,16 @@ namespace igl {
                     for (Index i = 0; i < cycleedges.size(); i++)
                     {
                         usededges.insert(cycleedges[i]);
-                    }                    
+                    }
                 }
             }
         }
     }
   }
 }
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template void igl::cut_to_disk<Eigen::Matrix<int, -1, -1, 0, -1, -1>, int>(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&);
+
+#endif

include/igl/tet_tet_adjacency.cpp

@@ -66,3 +66,8 @@ igl::tet_tet_adjacency(
   DerivedTT TTi;
   tet_tet_adjacency(T, TT, TTi);
 }
+
+#ifdef IGL_STATIC_LIBRARY
+template void igl::tet_tet_adjacency<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<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
+

tests/CMakeLists.txt

@@ -0,0 +1,81 @@
+cmake_minimum_required(VERSION 3.1)
+project(libigl_tutorials)
+message(STATUS "CMAKE_C_COMPILER: ${CMAKE_C_COMPILER}")
+message(STATUS "CMAKE_CXX_COMPILER: ${CMAKE_CXX_COMPILER}")
+
+### conditionally compile certain modules depending on libraries found on the system
+list(APPEND CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/../cmake)
+#find_package(CGAL QUIET COMPONENTS Core)
+find_package(MATLAB QUIET)
+find_package(MOSEK)
+
+### libIGL options: choose between header only and compiled static library
+option(LIBIGL_USE_STATIC_LIBRARY "Use LibIGL as static library" ON)
+option(LIBIGL_WITH_EMBREE      "Use Embree"         OFF)
+
+### libIGL options: choose your dependencies (by default everything is OFF, in this example we need the viewer)
+option(LIBIGL_USE_STATIC_LIBRARY     "Use libigl as static library" ON)
+option(LIBIGL_WITH_CGAL              "Use CGAL"           OFF)
+option(LIBIGL_WITH_COMISO            "Use CoMiso"         OFF)
+option(LIBIGL_WITH_CORK              "Use Cork"           OFF)
+option(LIBIGL_WITH_EMBREE            "Use Embree"         OFF)
+option(LIBIGL_WITH_MATLAB            "Use Matlab"         OFF)
+option(LIBIGL_WITH_MOSEK             "Use MOSEK"          OFF)
+option(LIBIGL_WITH_OPENGL            "Use OpenGL"         OFF)
+option(LIBIGL_WITH_OPENGL_GLFW       "Use GLFW"           OFF)
+option(LIBIGL_WITH_OPENGL_GLFW_IMGUI "Use ImGui"          OFF)
+option(LIBIGL_WITH_PNG               "Use PNG"            OFF)
+option(LIBIGL_WITH_TETGEN            "Use Tetgen"         OFF)
+option(LIBIGL_WITH_TRIANGLE          "Use Triangle"       OFF)
+option(LIBIGL_WITH_VIEWER            "Use OpenGL viewer"  OFF)
+option(LIBIGL_WITH_XML               "Use XML"            OFF)
+option(LIBIGL_WITH_PYTHON            "Use Python"         OFF)
+### End
+
+### Adding libIGL: choose the path to your local copy libIGL
+include(libigl)
+
+### Download data
+igl_download_test_data()
+set(IGL_TEST_DATA ${LIBIGL_EXTERNAL}/../tests/data)
+
+### Download Google unit test framework.
+igl_download_googletest()
+
+SET(TEST_ROOT_DIR ${CMAKE_CURRENT_LIST_DIR}) 
+INCLUDE_DIRECTORIES(${TEST_ROOT_DIR})
+
+# Set TEST_DIR definition
+ADD_DEFINITIONS(-DLIBIGL_DATA_DIR="${IGL_TEST_DATA}")
+
+# Add googletest googlemock support
+ENABLE_TESTING()
+ADD_SUBDIRECTORY(
+  ${LIBIGL_EXTERNAL}/googletest/googlemock
+  ${CMAKE_CURRENT_BINARY_DIR}/gtest)
+SET(GTEST_BOTH_LIBRARIES gtest gmock)
+INCLUDE_DIRECTORIES(${gmock_SOURCE_DIR})
+INCLUDE_DIRECTORIES(${gmock_SOURCE_DIR}/include)
+INCLUDE_DIRECTORIES(${gtest_SOURCE_DIR})
+INCLUDE_DIRECTORIES(${gtest_SOURCE_DIR}/include)
+
+# Process code in each subdirectories: add in decreasing order of complexity
+# (last added will run first and those should be the fastest tests)
+IF(LIBIGL_WITH_MOSEK)
+  ADD_SUBDIRECTORY(${TEST_ROOT_DIR}/include/igl/mosek)
+ENDIF()
+
+IF(LIBIGL_WITH_CGAL)
+  ADD_SUBDIRECTORY(${TEST_ROOT_DIR}/include/igl/copyleft/boolean)
+  ADD_SUBDIRECTORY(${TEST_ROOT_DIR}/include/igl/copyleft/cgal)
+ENDIF()
+
+IF(LIBIGL_WITH_TETGEN)
+  ADD_SUBDIRECTORY(${TEST_ROOT_DIR}/include/igl/copyleft/tetgen)
+ENDIF()
+
+IF(LIBIGL_WITH_COMISO)
+  ADD_SUBDIRECTORY(${TEST_ROOT_DIR}/include/igl/copyleft/comiso)
+ENDIF()
+
+ADD_SUBDIRECTORY(${TEST_ROOT_DIR}/include/igl)

tests/include/igl/CMakeLists.txt

@@ -0,0 +1,6 @@
+file(GLOB TEST_SRC_FILES *.cpp main.cpp)
+file(GLOB TEST_INC_FILES *.h *.inl)
+
+add_executable(igl_tests ${TEST_SRC_FILES} ${TEST_INC_FILES})
+target_link_libraries(igl_tests igl::core gtest_main)
+add_test(NAME run_igl_tests COMMAND igl_tests)

tests/include/igl/avg_edge_length.cpp

@@ -0,0 +1,62 @@
+#include <test_common.h>
+#include <igl/avg_edge_length.h>
+#include <iostream>
+
+
+TEST(avg_edge_length, cube)
+{
+  //The allowed error for this test
+  const double epsilon = 1e-15;
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  //This is a cube of dimensions 1.0x1.0x1.0
+  test_common::load_mesh("cube.obj", V, F);
+  //Create scaled versions of the cube
+  double scale = 1.0;
+  double huge_scale = 1.0e8;
+  double tiny_scale = 1.0e-8;
+  Eigen::MatrixXd V_huge = V * huge_scale;
+  Eigen::MatrixXd V_tiny = V * tiny_scale;
+  //Prepare another mesh with triangles along side diagonals of the cube
+  //These triangles are form a regular tetrahedron of side sqrt(2)
+  Eigen::MatrixXi F_tet(4,3);
+  F_tet << 4,6,1,
+            6,4,3,
+            4,1,3,
+            1,6,3;
+
+  //1. Check avg_edge_length function
+  double side_sq = 1.0; //squared lenght of a side
+  double diag_sq = 2.0; //squared lenght of a diagonal
+  double avg;
+
+  avg = igl::avg_edge_length(V,F);
+  ASSERT_NEAR((12.*sqrt(side_sq) + 6.*sqrt(diag_sq))/(12.+6.), avg, epsilon);
+
+  //Check the regular tetrahedron
+  avg = igl::avg_edge_length(V,F_tet);
+  ASSERT_NEAR(sqrt(diag_sq), avg, epsilon);
+
+
+  //Scale the cube to have huge sides
+  side_sq = huge_scale * huge_scale;  //squared lenght of a side
+  diag_sq = 2.0 * side_sq;  //squared lenght of a diagonal
+  avg = igl::avg_edge_length(V_huge,F);
+  ASSERT_NEAR((12.*sqrt(side_sq) + 6.*sqrt(diag_sq))/(12.+6.), avg, epsilon*huge_scale);
+
+  //Check the equilateral triangles
+  avg = igl::avg_edge_length(V_huge,F_tet);
+  ASSERT_NEAR(sqrt(diag_sq), avg, epsilon*huge_scale);
+
+
+  //Scale the cube to have tiny sides
+  side_sq = tiny_scale * tiny_scale;  //squared lenght of a side
+  diag_sq = 2.0 * side_sq;  //squared lenght of a diagonal
+  avg = igl::avg_edge_length(V_tiny,F);
+  ASSERT_NEAR((12.*sqrt(side_sq) + 6.*sqrt(diag_sq))/(12.+6.), avg, epsilon*tiny_scale);
+
+  //Check the regular tetrahedron
+  avg = igl::avg_edge_length(V_tiny,F_tet);
+  ASSERT_NEAR(sqrt(diag_sq), avg, epsilon*tiny_scale);
+
+}

tests/include/igl/bbw.cpp

@@ -0,0 +1,27 @@
+#include <gtest/gtest.h>
+#include <test_common.h>
+#include <igl/boundary_conditions.h>
+#include <igl/readMESH.h>
+#include <igl/writeDMAT.h>
+#include <igl/readTGF.h>
+#include <igl/bbw.h>
+
+TEST(bbw, decimated_knight)
+{
+  Eigen::MatrixXd V,C;
+  Eigen::MatrixXi T,F,E;
+  igl::readMESH(test_common::data_path("decimated-knight.mesh"),V,T,F);
+  igl::readTGF(test_common::data_path("decimated-knight.tgf"),C,E);
+  Eigen::MatrixXd W_groundtruth, Was, Wmo;
+  igl::readDMAT(
+    test_common::data_path("decimated-knight-matlab-active-set.dmat"),W_groundtruth);
+  Eigen::VectorXi b;
+  Eigen::MatrixXd bc;
+  igl::boundary_conditions(V,T,C,Eigen::VectorXi(),E,Eigen::MatrixXi(),b,bc);
+  igl::BBWData params;
+  params.active_set_params.max_iter = 100;
+  igl::bbw(V,T,b,bc,params,Was);
+  // igl::writeDMAT("decimated-knight-as.dmat",Was);
+  ASSERT_LT( (Was-W_groundtruth).array().abs().maxCoeff() ,1e-4);
+}
+

tests/include/igl/boundary_loop.cpp

@@ -0,0 +1,53 @@
+#include <test_common.h>
+#include <igl/boundary_loop.h>
+#include <vector>
+#include <algorithm>
+#include <iostream>
+
+TEST(boundary_loop, cube)
+{
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  //This is a cube of dimensions 1.0x1.0x1.0
+  test_common::load_mesh("cube.off", V, F);
+
+  //Compute Boundary Loop
+  Eigen::VectorXi boundary;
+  igl::boundary_loop(F, boundary);
+
+  //The cube has no boundary
+  ASSERT_EQ(0, boundary.size());
+}
+
+TEST(boundary_loop, bunny)
+{
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  //Load the Stanford bunny
+  test_common::load_mesh("bunny.off", V, F);
+
+  //Compute list of ordered boundary loops for a manifold mesh
+  std::vector<std::vector<int> >boundaries;
+  igl::boundary_loop(F, boundaries);
+
+  //Compare our result with known results taken from meshlab
+  ASSERT_EQ(5, boundaries.size());
+
+  //Compute min, max and sum of boundaries
+  size_t boundaryMin=9999999;
+  size_t boundaryMax=0;
+  size_t boundarySum=0;
+  for(size_t i=0; i<boundaries.size(); i++)
+  {
+      boundarySum += boundaries[i].size();
+      boundaryMax = std::max(boundaryMax, boundaries[i].size());
+      boundaryMin = std::min(boundaryMin, boundaries[i].size());
+  }
+
+  //Total boundary has 223 vertex
+  ASSERT_EQ(223, boundarySum);
+  //Largest loop has 80 vertex
+  ASSERT_EQ(80, boundaryMax);
+  //Smallest loop has 22 vertex
+  ASSERT_EQ(22, boundaryMin);
+}

tests/include/igl/copyleft/boolean/CMakeLists.txt

@@ -0,0 +1,7 @@
+# Check if CGAL is available
+file(GLOB TEST_SRC_FILES *.cpp main.cpp)
+file(GLOB TEST_INC_FILES *.h *.inl)
+
+add_executable(igl_boolean_tests ${TEST_SRC_FILES} ${TEST_INC_FILES})
+target_link_libraries(igl_boolean_tests igl::core igl::cgal gtest_main)
+add_test(NAME run_igl_boolean_tests COMMAND igl_boolean_tests)

tests/include/igl/copyleft/boolean/main.cpp

@@ -0,0 +1,6 @@
+#include <gtest/gtest.h>
+
+int main(int argc, char **argv) {
+      ::testing::InitGoogleTest(&argc, argv);
+        return RUN_ALL_TESTS();
+}

tests/include/igl/copyleft/boolean/mesh_boolean.cpp

@@ -0,0 +1,103 @@
+#include <test_common.h>
+
+#include <vector>
+
+#include <igl/copyleft/cgal/mesh_boolean.h>
+#include <igl/MeshBooleanType.h>
+#include <igl/exterior_edges.h>
+#include <igl/is_vertex_manifold.h>
+#include <igl/unique_edge_map.h>
+#include <igl/is_edge_manifold.h>
+
+namespace mesh_boolean_test {
+
+    template<typename DerivedF>
+    void assert_no_exterior_edges(
+            const Eigen::PlainObjectBase<DerivedF>& F) {
+        Eigen::MatrixXi Eb;
+        igl::exterior_edges(F, Eb);
+
+        ASSERT_EQ(0, Eb.rows());
+    }
+
+    template<typename DerivedV, typename DerivedF>
+    void assert_is_manifold(
+            const Eigen::PlainObjectBase<DerivedV>& V,
+            const Eigen::PlainObjectBase<DerivedF>& F) {
+        Eigen::MatrixXi B;
+        ASSERT_TRUE(igl::is_vertex_manifold(F, B));
+        ASSERT_TRUE(igl::is_edge_manifold(F));
+    }
+
+    template<typename DerivedV, typename DerivedF>
+    void assert_genus_eq(
+            const Eigen::PlainObjectBase<DerivedV>& V,
+            const Eigen::PlainObjectBase<DerivedF>& F,
+            const int genus) {
+        const int num_vertices = V.rows();
+        const int num_faces = F.rows();
+
+        Eigen::Matrix<
+            typename DerivedF::Scalar,
+            Eigen::Dynamic,
+            Eigen::Dynamic>
+            E, uE, EMAP;
+        std::vector<std::vector<size_t> > uE2E;
+        igl::unique_edge_map(F, E, uE, EMAP, uE2E);
+
+        const int num_edges = uE.rows();
+        const int euler = num_vertices - num_edges + num_faces;
+        ASSERT_EQ(euler, 2 - 2 * genus);
+    }
+
+TEST(MeshBoolean, TwoCubes) {
+    Eigen::MatrixXd V1;
+    Eigen::MatrixXi F1;
+    test_common::load_mesh("two-boxes-bad-self-union.ply", V1, F1);
+
+    Eigen::MatrixXd V2(0, 3);
+    Eigen::MatrixXi F2(0, 3);
+
+    Eigen::MatrixXd Vo;
+    Eigen::MatrixXi Fo;
+
+    igl::copyleft::cgal::mesh_boolean(V1, F1, V2, F2,
+            igl::MESH_BOOLEAN_TYPE_UNION,
+            Vo, Fo);
+
+    assert_no_exterior_edges(Fo);
+    assert_is_manifold(Vo, Fo);
+    assert_genus_eq(Vo, Fo, 0);
+}
+
+TEST(MeshBoolean, MinusTest) {
+    // Many thanks to Eric Yao for submitting this test case.
+    Eigen::MatrixXd V1, V2, Vo;
+    Eigen::MatrixXi F1, F2, Fo;
+    test_common::load_mesh("boolean_minus_test_cube.obj", V1, F1);
+    test_common::load_mesh("boolean_minus_test_green.obj", V2, F2);
+
+    igl::copyleft::cgal::mesh_boolean(V1, F1, V2, F2,
+            igl::MESH_BOOLEAN_TYPE_MINUS,
+            Vo, Fo);
+
+    assert_no_exterior_edges(Fo);
+    assert_is_manifold(Vo, Fo);
+    assert_genus_eq(Vo, Fo, 1);
+}
+
+TEST(MeshBoolean, IntersectWithSelf) {
+    Eigen::MatrixXd V1, Vo;
+    Eigen::MatrixXi F1, Fo;
+    test_common::load_mesh("cube.obj", V1, F1);
+
+    igl::copyleft::cgal::mesh_boolean(V1, F1, V1, F1,
+            igl::MESH_BOOLEAN_TYPE_INTERSECT,
+            Vo, Fo);
+
+    assert_no_exterior_edges(Fo);
+    assert_is_manifold(Vo, Fo);
+    assert_genus_eq(Vo, Fo, 0);
+}
+
+}

tests/include/igl/copyleft/cgal/CMakeLists.txt

@@ -0,0 +1,6 @@
+file(GLOB TEST_SRC_FILES *.cpp main.cpp)
+file(GLOB TEST_INC_FILES *.h *.inl)
+
+add_executable(igl_cgal_tests ${TEST_SRC_FILES} ${TEST_INC_FILES})
+target_link_libraries(igl_cgal_tests igl::core igl::cgal gtest_main)
+add_test(NAME run_igl_cgal_tests COMMAND igl_cgal_tests)

tests/include/igl/copyleft/cgal/CSGTree.cpp

@@ -0,0 +1,17 @@
+#include <test_common.h>
+
+#include <igl/copyleft/cgal/CSGTree.h>
+
+TEST(CSGTree, extrusion) {
+    Eigen::MatrixXd V;
+    Eigen::MatrixXi F;
+    test_common::load_mesh("extrusion.obj", V, F);
+    igl::copyleft::cgal::CSGTree tree(V, F);
+    igl::copyleft::cgal::CSGTree inter(tree, tree, "i"); // returns error
+
+    Eigen::MatrixXd V2 = inter.cast_V<Eigen::MatrixXd>();
+    Eigen::MatrixXi F2 = inter.F();
+
+    ASSERT_EQ(V.rows(), V2.rows());
+    ASSERT_EQ(F.rows(), F2.rows());
+}

tests/include/igl/copyleft/cgal/hausdorff.cpp

@@ -0,0 +1,55 @@
+#include <test_common.h>
+
+#include <CGAL/Simple_cartesian.h>
+#include <igl/copyleft/cgal/hausdorff.h>
+#include <igl/copyleft/cgal/point_mesh_squared_distance.h>
+#include <igl/upsample.h>
+
+TEST(hausdorff, knightVScheburashka) 
+{
+  Eigen::MatrixXd VA,VB;
+  Eigen::MatrixXi FA,FB;
+  test_common::load_mesh("decimated-knight.obj", VA, FA);
+  test_common::load_mesh("cheburashka.off", VB, FB);
+  //typedef CGAL::Epeck Kernel;
+  typedef CGAL::Simple_cartesian<double> Kernel;
+  CGAL::AABB_tree<
+    CGAL::AABB_traits<Kernel, 
+      CGAL::AABB_triangle_primitive<Kernel, 
+        typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
+      >
+    >
+  > treeB;
+  std::vector<CGAL::Triangle_3<Kernel> > TB;
+  {
+    igl::copyleft::cgal::point_mesh_squared_distance_precompute(VB,FB,treeB,TB);
+  }
+  std::vector<Eigen::VectorXd> U(5);
+  for(int j = 0;j<U.size();j++)
+  {
+    if(j>0)
+    {
+      igl::upsample(Eigen::MatrixXd(VA),Eigen::MatrixXi(FA),VA,FA);
+    }
+    const int m = FA.rows();
+    U[j].resize(m);
+    for(int i = 0;i<m;i++)
+    {
+      Eigen::MatrixXd Vi(3,3);
+      Vi<<VA.row(FA(i,0)),VA.row(FA(i,1)),VA.row(FA(i,2));
+      double l;
+      igl::copyleft::cgal::hausdorff(Vi,treeB,TB,l,U[j](i));
+      if(j>0&&i%4==3)
+      {
+        double u4 = std::numeric_limits<double>::infinity();
+        for(int u = 0;u<4;u++)
+        {
+          u4 = std::min(u4,U[j](i-u));
+        }
+        ASSERT_LE(u4,U[j-1](i/4));
+      }
+    }
+    break;
+  }
+}
+

tests/include/igl/copyleft/cgal/main.cpp

@@ -0,0 +1,6 @@
+#include <gtest/gtest.h>
+
+int main(int argc, char **argv) {
+      ::testing::InitGoogleTest(&argc, argv);
+        return RUN_ALL_TESTS();
+}

tests/include/igl/copyleft/cgal/order_facets_around_edges.cpp

@@ -0,0 +1,189 @@
+#include <test_common.h>
+
+#include <algorithm>
+#include <iostream>
+#include <vector>
+
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <igl/copyleft/cgal/order_facets_around_edges.h>
+#include <igl/unique_edge_map.h>
+#include <igl/readDMAT.h>
+#include <igl/per_face_normals.h>
+
+namespace order_facets_around_edges_test {
+
+typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
+
+template<typename T>
+size_t index_of(const std::vector<T>& array, T val) {
+    auto loc = std::find(array.begin(), array.end(), val);
+    assert(loc != array.end());
+    return loc - array.begin();
+}
+
+void assert_consistently_oriented(size_t num_faces,
+        const std::vector<int>& expected_face_order,
+        const std::vector<int>& e_order) {
+    const size_t num_items = expected_face_order.size();
+    ASSERT_EQ(num_items, e_order.size());
+
+    std::vector<int> order(num_items);
+    std::transform(e_order.begin(), e_order.end(), order.begin(),
+            [=](int val) { return val % num_faces; });
+
+    size_t ref_start = index_of(order, expected_face_order[0]);
+    for (size_t i=0; i<num_items; i++) {
+        ASSERT_EQ(expected_face_order[i], order[(ref_start + i) % num_items]);
+    }
+}
+
+template<typename DerivedV, typename DerivedF>
+void assert_order(
+        const Eigen::PlainObjectBase<DerivedV>& V,
+        const Eigen::PlainObjectBase<DerivedF>& F,
+        size_t v0, size_t v1,
+        std::vector<int> expected_order, const std::string& normal="") {
+    Eigen::MatrixXi E, uE, EMAP;
+    std::vector<std::vector<int> > uE2E;
+    igl::unique_edge_map(F, E, uE, EMAP, uE2E);
+
+    std::vector<std::vector<int> > uE2oE;
+    std::vector<std::vector<bool> > uE2C;
+
+    if (normal != "") {
+        Eigen::MatrixXd N;
+        //igl::per_face_normals_stable(V, F, N);
+        //igl::per_face_normals(V, F, N);
+        test_common::load_matrix(normal, N);
+        igl::copyleft::cgal::order_facets_around_edges(
+          V, F, N, uE, uE2E, uE2oE, uE2C);
+    } else {
+        igl::copyleft::cgal::order_facets_around_edges(
+          V, F, uE, uE2E, uE2oE, uE2C);
+    }
+
+    const size_t num_uE = uE.rows();
+    for (size_t i=0; i<num_uE; i++) {
+        const auto& order = uE2oE[i];
+        const auto& cons  = uE2C[i];
+        Eigen::VectorXi e = uE.row(i);
+        if (order.size() <= 1) continue;
+        if (e[0] != v0 && e[0] != v1) continue;
+        if (e[1] != v0 && e[1] != v1) continue;
+        if (e[0] == v1 && e[1] == v0) {
+            std::reverse(expected_order.begin(), expected_order.end());
+        }
+        assert_consistently_oriented(F.rows(), expected_order, order);
+    }
+}
+
+TEST(copyleft_cgal_order_facets_around_edges, Simple) {
+    Eigen::MatrixXd V(4, 3);
+    V << 0.0, 0.0, 0.0,
+         1.0, 0.0, 0.0,
+         0.0, 1.0, 0.0,
+         1.0, 1.0, 0.0;
+    Eigen::MatrixXi F(2, 3);
+    F << 0, 1, 2,
+         2, 1, 3;
+
+    assert_order(V, F, 1, 2, {0, 1});
+}
+
+TEST(copyleft_cgal_order_facets_around_edges, TripletFaces) {
+    Eigen::MatrixXd V(5, 3);
+    V << 0.0, 0.0, 0.0,
+         1.0, 0.0, 0.0,
+         0.0, 1.0, 0.0,
+         1.0, 1.0, 0.0,
+         0.0, 0.0, 1.0;
+    Eigen::MatrixXi F(3, 3);
+    F << 0, 1, 2,
+         2, 1, 3,
+         1, 2, 4;
+
+    assert_order(V, F, 1, 2, {0, 1, 2});
+}
+
+TEST(copyleft_cgal_order_facets_around_edges, DuplicatedFaces) {
+    Eigen::MatrixXd V(5, 3);
+    V << 0.0, 0.0, 0.0,
+         1.0, 0.0, 0.0,
+         0.0, 1.0, 0.0,
+         1.0, 1.0, 0.0,
+         0.0, 0.0, 1.0;
+    Eigen::MatrixXi F(4, 3);
+    F << 0, 1, 2,
+         2, 1, 3,
+         1, 2, 4,
+         4, 1, 2;
+
+    assert_order(V, F, 1, 2, {0, 1, 3, 2});
+}
+
+TEST(copyleft_cgal_order_facets_around_edges, MultipleDuplicatedFaces) {
+    Eigen::MatrixXd V(5, 3);
+    V << 0.0, 0.0, 0.0,
+         1.0, 0.0, 0.0,
+         0.0, 1.0, 0.0,
+         1.0, 1.0, 0.0,
+         0.0, 0.0, 1.0;
+    Eigen::MatrixXi F(6, 3);
+    F << 0, 1, 2,
+         1, 2, 0,
+         2, 1, 3,
+         1, 3, 2,
+         1, 2, 4,
+         4, 1, 2;
+
+    assert_order(V, F, 1, 2, {1, 0, 2, 3, 5, 4});
+}
+
+TEST(copyleft_cgal_order_facets_around_edges, Debug) {
+    Eigen::MatrixXd V(5, 3);
+    V <<
+        -44.3205080756887781, 4.22994972382184579e-15, 75,
+        -27.933756729740665, -48.382685902179837, 75,
+        -55.8675134594812945, -2.81996648254789745e-15, 75,
+        -27.933756729740665, -48.382685902179837, 70,
+        -31.4903810567666049, -42.2224318643354408, 85;
+
+    Eigen::MatrixXi F(3, 3);
+    F << 1, 0, 2,
+         2, 3, 1,
+         4, 1, 2;
+
+    assert_order(V, F, 1, 2, {0, 2, 1});
+}
+
+TEST(copyleft_cgal_order_facets_around_edges, Debug2) {
+    Eigen::MatrixXd V(5, 3);
+    V <<
+        -22.160254037844382, 38.3826859021798441, 75,
+        -27.9337567297406331, 48.3826859021798654, 75,
+        27.9337567297406544, 48.3826859021798512, 75,
+        27.9337567297406544, 48.3826859021798512, 70,
+        20.8205080756887924, 48.3826859021798512, 85;
+    Eigen::MatrixXi F(3, 3);
+    F << 1, 0, 2,
+         3, 1, 2,
+         2, 4, 1;
+
+    assert_order(V, F, 1, 2, {1, 0, 2});
+}
+
+TEST(copyleft_cgal_order_facets_around_edges, NormalSensitivity) {
+    // This example shows that epsilon difference in normal vectors could
+    // results in very different ordering of facets.
+
+    Eigen::MatrixXd V;
+    test_common::load_matrix("duplicated_faces_V.dmat", V);
+    Eigen::MatrixXi F;
+    test_common::load_matrix("duplicated_faces_F.dmat", F);
+
+    assert_order(V, F, 223, 224, {2, 0, 3, 1}, "duplicated_faces_N1.dmat");
+    assert_order(V, F, 223, 224, {0, 3, 2, 1}, "duplicated_faces_N2.dmat");
+}
+
+
+}

tests/include/igl/copyleft/cgal/outer_facet.cpp

@@ -0,0 +1,108 @@
+#include <test_common.h>
+
+#include <igl/copyleft/cgal/outer_facet.h>
+
+namespace OuterFacetHelper {
+
+/**
+ * Check if the outer facet is indeed valid.
+ * Assumption: mesh is closed.
+ */
+template<typename DerivedV, typename DerivedF>
+void assert_outer_facet_is_correct(
+        const Eigen::PlainObjectBase<DerivedV>& V,
+        const Eigen::PlainObjectBase<DerivedF>& F,
+        size_t fid, bool flipped) {
+    // Todo.
+}
+
+TEST(OuterFacet, Simple) {
+    Eigen::MatrixXd V;
+    Eigen::MatrixXi F;
+    test_common::load_mesh("cube.obj", V, F);
+
+    const size_t num_faces = F.rows();
+
+    Eigen::VectorXi I(num_faces);
+    I.setLinSpaced(num_faces, 0, num_faces-1);
+
+    size_t fid = num_faces + 1;
+    bool flipped;
+    igl::copyleft::cgal::outer_facet(V, F, I, fid, flipped);
+
+    ASSERT_LT(fid, num_faces);
+    ASSERT_FALSE(flipped);
+}
+
+TEST(OuterFacet, DuplicatedOppositeFaces) {
+    Eigen::MatrixXd V;
+    Eigen::MatrixXi F1;
+    test_common::load_mesh("cube.obj", V, F1);
+
+    Eigen::MatrixXi F2 = F1;
+    F2.col(0).swap(F2.col(1));
+
+    Eigen::MatrixXi F(F1.rows()*2, F1.cols());
+    F << F1, F2;
+
+    Eigen::VectorXi I(F.rows());
+    I.setLinSpaced(F.rows(), 0, F.rows()-1);
+
+    size_t fid = F.rows() + 1;
+    bool flipped;
+    igl::copyleft::cgal::outer_facet(V, F, I, fid, flipped);
+
+    ASSERT_LT(fid, F.rows());
+    ASSERT_FALSE(flipped);
+}
+
+TEST(OuterFacet, FullyDegnerated) {
+    Eigen::MatrixXd V;
+    Eigen::MatrixXi F;
+    test_common::load_mesh("degenerated.obj", V, F);
+
+    Eigen::VectorXi I(F.rows());
+    I.setLinSpaced(F.rows(), 0, F.rows()-1);
+
+    size_t fid = F.rows() + 1;
+    bool flipped;
+    igl::copyleft::cgal::outer_facet(V, F, I, fid, flipped);
+
+    ASSERT_LT(fid, F.rows());
+    ASSERT_FALSE(flipped);
+}
+
+TEST(OuterFacet, InvertedNormal) {
+    Eigen::MatrixXd V;
+    Eigen::MatrixXi F;
+    test_common::load_mesh("cube.obj", V, F);
+    F.col(0).swap(F.col(1));
+
+    Eigen::VectorXi I(F.rows());
+    I.setLinSpaced(F.rows(), 0, F.rows()-1);
+
+    size_t fid = F.rows() + 1;
+    bool flipped;
+    igl::copyleft::cgal::outer_facet(V, F, I, fid, flipped);
+
+    ASSERT_LT(fid, F.rows());
+    ASSERT_TRUE(flipped);
+}
+
+TEST(OuterFacet, SliverTet) {
+    Eigen::MatrixXd V;
+    Eigen::MatrixXi F;
+    test_common::load_mesh("sliver_tet.ply", V, F);
+
+    Eigen::VectorXi I(F.rows());
+    I.setLinSpaced(F.rows(), 0, F.rows()-1);
+
+    size_t fid = F.rows() + 1;
+    bool flipped;
+    igl::copyleft::cgal::outer_facet(V, F, I, fid, flipped);
+
+    ASSERT_LT(fid, F.rows());
+    ASSERT_FALSE(flipped);
+}
+
+}

tests/include/igl/copyleft/cgal/outer_hull.cpp

@@ -0,0 +1,14 @@
+#include <test_common.h>
+
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <igl/copyleft/cgal/outer_hull.h>
+
+TEST(OuterHull, CubeWithFold) {
+    Eigen::MatrixXd V;
+    Eigen::MatrixXi F;
+    test_common::load_mesh("cube_with_fold.ply", V, F);
+
+    Eigen::MatrixXi G,J,flip;
+    // Is this just checking that it doesn't crash?
+    igl::copyleft::cgal::outer_hull_legacy(V, F, G, J, flip);
+}

tests/include/igl/copyleft/cgal/peel_outer_hull_layers.cpp

@@ -0,0 +1,78 @@
+#include <test_common.h>
+#include <iostream>
+#include <Eigen/Dense>
+
+#include <igl/copyleft/cgal/peel_outer_hull_layers.h>
+#include <igl/copyleft/cgal/remesh_self_intersections.h>
+#include <igl/copyleft/cgal/RemeshSelfIntersectionsParam.h>
+#include <igl/per_face_normals.h>
+#include <igl/remove_unreferenced.h>
+#include <igl/writeOBJ.h>
+
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+
+TEST(copyleft_cgal_peel_outer_hull_layers, TwoCubes) {
+    Eigen::MatrixXd V;
+    Eigen::MatrixXi F;
+    test_common::load_mesh("two-boxes-bad-self-union.ply", V, F);
+    ASSERT_EQ(486, V.rows());
+    ASSERT_EQ(708, F.rows());
+
+    typedef CGAL::Exact_predicates_exact_constructions_kernel K;
+    typedef K::FT Scalar;
+    typedef Eigen::Matrix<Scalar,
+            Eigen::Dynamic,
+            Eigen::Dynamic> MatrixXe;
+
+    MatrixXe Vs;
+    Eigen::MatrixXi Fs, IF;
+    Eigen::VectorXi J, IM;
+    igl::copyleft::cgal::RemeshSelfIntersectionsParam param;
+    igl::copyleft::cgal::remesh_self_intersections(V, F, param, Vs, Fs, IF, J, IM);
+
+    std::for_each(Fs.data(),Fs.data()+Fs.size(),
+            [&IM](int & a){ a=IM(a); });
+    MatrixXe Vt;
+    Eigen::MatrixXi Ft;
+    igl::remove_unreferenced(Vs,Fs,Vt,Ft,IM);
+    const size_t num_faces = Ft.rows();
+
+    Eigen::VectorXi I, flipped;
+    size_t num_peels = igl::copyleft::cgal::peel_outer_hull_layers(Vt, Ft, I, flipped);
+
+    Eigen::MatrixXd vertices(Vt.rows(), Vt.cols());
+    std::transform(Vt.data(), Vt.data() + Vt.rows() * Vt.cols(),
+            vertices.data(), [](Scalar v) { return CGAL::to_double(v); });
+    igl::writeOBJ("debug.obj", vertices, Ft);
+
+    ASSERT_EQ(num_faces, I.rows());
+    ASSERT_EQ(0, I.minCoeff());
+    ASSERT_EQ(1, I.maxCoeff());
+}
+
+TEST(PeelOuterHullLayers, CubeWithFold) {
+    Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> V;
+    Eigen::MatrixXi F;
+    test_common::load_mesh("cube_with_fold.ply", V, F);
+
+    typedef CGAL::Exact_predicates_exact_constructions_kernel K;
+    typedef K::FT Scalar;
+    typedef Eigen::Matrix<Scalar,
+            Eigen::Dynamic,
+            Eigen::Dynamic> MatrixXe;
+
+    MatrixXe Vs;
+    Eigen::MatrixXi Fs, IF;
+    Eigen::VectorXi J, IM;
+    igl::copyleft::cgal::RemeshSelfIntersectionsParam param;
+    igl::copyleft::cgal::remesh_self_intersections(V, F, param, Vs, Fs, IF, J, IM);
+
+    std::for_each(Fs.data(),Fs.data()+Fs.size(),
+            [&IM](int & a){ a=IM(a); });
+    MatrixXe Vt;
+    Eigen::MatrixXi Ft;
+    igl::remove_unreferenced(Vs,Fs,Vt,Ft,IM);
+
+    Eigen::VectorXi I, flipped;
+    size_t num_peels = igl::copyleft::cgal::peel_outer_hull_layers(Vt, Ft, I, flipped);
+}

tests/include/igl/copyleft/cgal/points_inside_component.cpp

@@ -0,0 +1,84 @@
+#include <test_common.h>
+
+#include <igl/copyleft/cgal/points_inside_component.h>
+#include <limits>
+
+namespace PointInsideComponentHelper {
+
+TEST(PointInsideComponent, simple) {
+    Eigen::MatrixXd V1;
+    Eigen::MatrixXi F1;
+    test_common::load_mesh("cube.obj", V1, F1);
+
+    Eigen::MatrixXd P(4, 3);
+    P << 0.0, 0.0, 0.0,
+         1.0, 0.0, 0.0,
+         0.0, 1.0, 0.0,
+         0.0, 0.0, 1.0;
+    Eigen::VectorXi inside;
+
+    EXPECT_NO_THROW(igl::copyleft::cgal::points_inside_component(V1, F1, P, inside));
+    ASSERT_EQ(1, inside[0]);
+    ASSERT_EQ(0, inside[1]);
+    ASSERT_EQ(0, inside[2]);
+    ASSERT_EQ(0, inside[3]);
+}
+
+TEST(PointInsideComponent, near_boundary) {
+    Eigen::MatrixXd V1;
+    Eigen::MatrixXi F1;
+    test_common::load_mesh("cube.obj", V1, F1);
+
+    const double EPS = std::numeric_limits<double>::epsilon();
+    Eigen::MatrixXd P(6, 3);
+    P << 0.5 + EPS, 0.0, 0.0,
+         0.0, 0.5 + EPS, 0.0,
+         0.0, 0.0, 0.5 + EPS,
+         0.5 - EPS, 0.0, 0.0,
+         0.0, 0.5 - EPS, 0.0,
+         0.0, 0.0, 0.5 - EPS;
+
+    Eigen::VectorXi inside;
+    EXPECT_NO_THROW(igl::copyleft::cgal::points_inside_component(V1, F1, P, inside));
+    ASSERT_EQ(0, inside[0]);
+    ASSERT_EQ(0, inside[1]);
+    ASSERT_EQ(0, inside[2]);
+    ASSERT_EQ(1, inside[3]);
+    ASSERT_EQ(1, inside[4]);
+    ASSERT_EQ(1, inside[5]);
+}
+
+TEST(PointInsideComponent, near_corner) {
+    Eigen::MatrixXd V1;
+    Eigen::MatrixXi F1;
+    test_common::load_mesh("cube.obj", V1, F1);
+
+    const double EPS = std::numeric_limits<double>::epsilon();
+    Eigen::MatrixXd P_out(8, 3);
+    P_out << 0.5 + EPS, 0.5 + EPS, 0.5 + EPS,
+            -0.5 - EPS, 0.5 + EPS, 0.5 + EPS,
+             0.5 + EPS,-0.5 - EPS, 0.5 + EPS,
+            -0.5 - EPS,-0.5 - EPS, 0.5 + EPS,
+             0.5 + EPS, 0.5 + EPS,-0.5 - EPS,
+            -0.5 - EPS, 0.5 + EPS,-0.5 - EPS,
+             0.5 + EPS,-0.5 - EPS,-0.5 - EPS,
+            -0.5 - EPS,-0.5 - EPS,-0.5 - EPS;
+
+    Eigen::VectorXi inside;
+    EXPECT_NO_THROW(igl::copyleft::cgal::points_inside_component(V1, F1, P_out, inside));
+    ASSERT_TRUE((inside.array()==0).all());
+
+    Eigen::MatrixXd P_in(8, 3);
+    P_in << 0.5 - EPS, 0.5 - EPS, 0.5 - EPS,
+           -0.5 + EPS, 0.5 - EPS, 0.5 - EPS,
+            0.5 - EPS,-0.5 + EPS, 0.5 - EPS,
+           -0.5 + EPS,-0.5 + EPS, 0.5 - EPS,
+            0.5 - EPS, 0.5 - EPS,-0.5 + EPS,
+           -0.5 + EPS, 0.5 - EPS,-0.5 + EPS,
+            0.5 - EPS,-0.5 + EPS,-0.5 + EPS,
+           -0.5 + EPS,-0.5 + EPS,-0.5 + EPS;
+    EXPECT_NO_THROW(igl::copyleft::cgal::points_inside_component(V1, F1, P_in, inside));
+    ASSERT_TRUE((inside.array()==1).all());
+}
+
+}

tests/include/igl/copyleft/cgal/remesh_self_intersections.cpp

@@ -0,0 +1,22 @@
+#include <test_common.h>
+#include <Eigen/Dense>
+
+#include <igl/copyleft/cgal/remesh_self_intersections.h>
+#include <igl/copyleft/cgal/RemeshSelfIntersectionsParam.h>
+
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+
+TEST(RemeshSelfIntersections, CubeWithFold) {
+    Eigen::MatrixXd V;
+    Eigen::MatrixXi F;
+    test_common::load_mesh("cube_with_fold.ply", V, F);
+
+    typedef CGAL::Exact_predicates_exact_constructions_kernel K;
+    typedef Eigen::Matrix<K::FT, Eigen::Dynamic, Eigen::Dynamic> MatrixXe;
+
+    MatrixXe VV;
+    Eigen::MatrixXi FF, IF;
+    Eigen::VectorXi J, IM;
+    igl::copyleft::cgal::RemeshSelfIntersectionsParam param;
+    igl::copyleft::cgal::remesh_self_intersections(V, F, param, VV, FF, IF, J, IM);
+}

tests/include/igl/copyleft/comiso/CMakeLists.txt

@@ -0,0 +1,6 @@
+file(GLOB TEST_SRC_FILES *.cpp main.cpp)
+file(GLOB TEST_INC_FILES *.h *.inl)
+
+add_executable(igl_comiso_tests ${TEST_SRC_FILES} ${TEST_INC_FILES})
+target_link_libraries(igl_comiso_tests igl::core igl::comiso gtest_main)
+add_test(NAME run_igl_comiso_tests COMMAND igl_comiso_tests)

tests/include/igl/copyleft/comiso/main.cpp

@@ -0,0 +1,6 @@
+#include <gtest/gtest.h>
+
+int main(int argc, char **argv) {
+      ::testing::InitGoogleTest(&argc, argv);
+        return RUN_ALL_TESTS();
+}

tests/include/igl/copyleft/comiso/miq.cpp

@@ -0,0 +1,131 @@
+#include <gtest/gtest.h>
+#include <test_common.h>
+#include <igl/avg_edge_length.h>
+#include <igl/barycenter.h>
+#include <igl/comb_cross_field.h>
+#include <igl/comb_frame_field.h>
+#include <igl/compute_frame_field_bisectors.h>
+#include <igl/cross_field_missmatch.h>
+#include <igl/cut_mesh_from_singularities.h>
+#include <igl/find_cross_field_singularities.h>
+#include <igl/local_basis.h>
+#include <igl/readOFF.h>
+#include <igl/rotate_vectors.h>
+#include <igl/copyleft/comiso/miq.h>
+#include <igl/copyleft/comiso/nrosy.h>
+#include <igl/PI.h>
+#include <igl/serialize.h>
+#include <sstream>
+#include <igl/writeDMAT.h>
+
+TEST(miq, 3_holes)
+{
+using namespace Eigen;
+
+// Input mesh
+Eigen::MatrixXd V;
+Eigen::MatrixXi F;
+
+// Face barycenters
+Eigen::MatrixXd B;
+
+// Cross field
+Eigen::MatrixXd X1,X2;
+
+// Bisector field
+Eigen::MatrixXd BIS1, BIS2;
+
+// Combed bisector
+Eigen::MatrixXd BIS1_combed, BIS2_combed;
+
+// Per-corner, integer mismatches
+Eigen::Matrix<int, Eigen::Dynamic, 3> MMatch;
+
+// Field singularities
+Eigen::Matrix<int, Eigen::Dynamic, 1> isSingularity, singularityIndex;
+
+// Per corner seams
+Eigen::Matrix<int, Eigen::Dynamic, 3> Seams;
+
+// Combed field
+Eigen::MatrixXd X1_combed, X2_combed;
+
+// Global parametrization
+Eigen::MatrixXd UV;
+Eigen::MatrixXi FUV;
+
+// Global parametrization (reference)
+Eigen::MatrixXd UV_ref;
+Eigen::MatrixXi FUV_ref;
+
+// Load a mesh in OFF format
+igl::readOFF(test_common::data_path("3holes.off"), V, F);
+
+double gradient_size = 50;
+double iter = 0;
+double stiffness = 5.0;
+bool direct_round = 0;
+
+// Compute face barycenters
+igl::barycenter(V, F, B);
+
+// Contrain one face
+VectorXi b(1);
+b << 0;
+MatrixXd bc(1, 3);
+bc << 1, 0, 0;
+
+// Create a smooth 4-RoSy field
+VectorXd S;
+igl::copyleft::comiso::nrosy(V, F, b, bc, VectorXi(), VectorXd(), MatrixXd(), 4, 0.5, X1, S);
+
+// Find the orthogonal vector
+MatrixXd B1, B2, B3;
+igl::local_basis(V, F, B1, B2, B3);
+X2 = igl::rotate_vectors(X1, VectorXd::Constant(1, igl::PI / 2), B1, B2);
+
+// Always work on the bisectors, it is more general
+igl::compute_frame_field_bisectors(V, F, X1, X2, BIS1, BIS2);
+
+// Comb the field, implicitly defining the seams
+igl::comb_cross_field(V, F, BIS1, BIS2, BIS1_combed, BIS2_combed);
+
+// Find the integer mismatches
+igl::cross_field_missmatch(V, F, BIS1_combed, BIS2_combed, true, MMatch);
+
+// Find the singularities
+igl::find_cross_field_singularities(V, F, MMatch, isSingularity, singularityIndex);
+
+// Cut the mesh, duplicating all vertices on the seams
+igl::cut_mesh_from_singularities(V, F, MMatch, Seams);
+
+// Comb the frame-field accordingly
+igl::comb_frame_field(V, F, X1, X2, BIS1_combed, BIS2_combed, X1_combed, X2_combed);
+
+// Global parametrization
+igl::copyleft::comiso::miq(V,
+          F,
+          X1_combed,
+          X2_combed,
+          MMatch,
+          isSingularity,
+          Seams,
+          UV,
+          FUV,
+          gradient_size,
+          stiffness,
+          direct_round,
+          iter,
+          5,
+          true);
+
+  // Refresh the test data
+  // igl::writeDMAT(test_common::data_path("3holes-miq-UV.dmat"),UV);
+  // igl::writeDMAT(test_common::data_path("3holes-miq-FUV.dmat"),FUV);
+
+  igl::readDMAT(test_common::data_path("3holes-miq-UV.dmat"),UV_ref);
+  igl::readDMAT(test_common::data_path("3holes-miq-FUV.dmat"),FUV_ref);
+
+  ASSERT_LT((UV-UV_ref).array().abs().maxCoeff() ,1e-6);
+  ASSERT_LT((FUV-FUV_ref).array().abs().maxCoeff() ,1e-6);
+}

tests/include/igl/copyleft/tetgen/CMakeLists.txt

@@ -0,0 +1,6 @@
+file(GLOB TEST_SRC_FILES *.cpp main.cpp)
+file(GLOB TEST_INC_FILES *.h *.inl)
+
+add_executable(igl_tetgen_tests ${TEST_SRC_FILES} ${TEST_INC_FILES})
+target_link_libraries(igl_tetgen_tests igl::core igl::tetgen gtest_main)
+add_test(NAME run_igl_tetgen_tests COMMAND igl_tetgen_tests)

tests/include/igl/copyleft/tetgen/main.cpp

@@ -0,0 +1,6 @@
+#include <gtest/gtest.h>
+
+int main(int argc, char **argv) {
+      ::testing::InitGoogleTest(&argc, argv);
+        return RUN_ALL_TESTS();
+}

tests/include/igl/copyleft/tetgen/tetrahedralize.cpp

@@ -0,0 +1,49 @@
+#include <test_common.h>
+
+#include <igl/setxor.h>
+#include <igl/copyleft/tetgen/tetrahedralize.h>
+
+TEST(tetrahedralize, single_tet) {
+  const Eigen::MatrixXd V = (Eigen::MatrixXd(4,3)<<
+    0,0,0,
+    1,0,0,
+    0,1,0,
+    0,0,1).finished();
+  Eigen::MatrixXi F(0,3);
+  Eigen::MatrixXd TV;
+  Eigen::MatrixXi TT,TF;
+  igl::copyleft::tetgen::tetrahedralize(V,F,"cpQ",TV,TT,TF);
+  ASSERT_EQ(TV.rows() , 4);
+  ASSERT_EQ(TT.rows() , 1);
+  ASSERT_EQ(TF.rows() , 4);
+  Eigen::MatrixXi TTcorrect = (Eigen::MatrixXi(1,4)<<0,1,2,3).finished();
+  {
+    Eigen::VectorXi TT_XOR,IA,IB;
+    igl::setxor(TT,TTcorrect,TT_XOR,IA,IB);
+    ASSERT_EQ(TT_XOR.size(),0);
+  }
+  test_common::assert_eq(TV,V);
+}
+
+TEST(tetrahedralize, schoenhardt) {
+  const Eigen::MatrixXd V = (Eigen::MatrixXd(6,3)<<
+    -246.2,-43.412,500,
+     160.7,-191.51,500,
+    85.505, 234.92,500,
+       250,      0,  0,
+      -125, 216.51,  0,
+      -125,-216.51,  0).finished();
+  const Eigen::MatrixXi F = (Eigen::MatrixXi(8,3)<<
+    1,2,0,
+    5,0,2,
+    3,1,0,
+    4,2,1,
+    0,5,3,
+    1,3,4,
+    2,4,5,
+    3,5,4).finished();
+  Eigen::MatrixXd TV;
+  Eigen::MatrixXi TT,TF;
+  igl::copyleft::tetgen::tetrahedralize(V,F,"pQ",TV,TT,TF);
+  ASSERT_GE(TV.rows() , V.rows());
+}

tests/include/igl/cotmatrix.cpp

@@ -0,0 +1,183 @@
+#include <test_common.h>
+#include <igl/cotmatrix.h>
+
+class cotmatrix : public ::testing::TestWithParam<std::string> {};
+
+TEST_P(cotmatrix, constant_in_null_space)
+{
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  Eigen::SparseMatrix<double> L;
+  // Load example mesh: GetParam() will be name of mesh file
+  test_common::load_mesh(GetParam(), V, F);
+  igl::cotmatrix(V,F,L);
+  ASSERT_EQ(V.rows(),L.rows());
+  ASSERT_EQ(L.rows(),L.cols());
+  Eigen::VectorXd C = Eigen::VectorXd::Ones(L.rows());
+  Eigen::VectorXd Z = Eigen::VectorXd::Zero(L.rows());
+  // ASSERT_EQ(a,b);
+  // ASSERT_TRUE(a==b);
+  // ASSERT_NEAR(a,b,1e-15)
+  ASSERT_LT(((L*C)-(Z)).norm(),1e-12);
+}
+
+INSTANTIATE_TEST_CASE_P
+(
+ all_meshes,
+ cotmatrix,
+ ::testing::ValuesIn(test_common::all_meshes()),
+ test_common::string_test_name
+);
+
+TEST(cotmatrix, cube)
+{
+  //The allowed error for this test
+  const double epsilon = 1e-15;
+
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  //This is a cube of dimensions 1.0x1.0x1.0
+  test_common::load_mesh("cube.obj", V, F);
+
+  //Scale the cube to have huge sides
+  Eigen::MatrixXd V_huge = V * 1.0e8;
+
+  //Scale the cube to have tiny sides
+  Eigen::MatrixXd V_tiny = V * 1.0e-8;
+
+  //Check cotmatrix (Laplacian)
+  //The laplacian for the cube is quite singular.
+  //Each edge in a diagonal has two opposite angles of 90, with cotangent 0.0 each
+  //Each edge in a side has two opposite angle of 45, with (half)cotangen 0.5 each
+  //So the cotangent matrix always are (0+0) or (0.5+0.5)
+  Eigen::SparseMatrix<double> L1;
+  igl::cotmatrix(V,F,L1);
+  ASSERT_EQ(V.rows(), L1.rows());
+  ASSERT_EQ(V.rows(), L1.cols());
+  for(int f = 0;f<L1.rows();f++)
+  {
+#ifdef IGL_EDGE_LENGTHS_SQUARED_H
+    //Hard assert if we have edge_lenght_squared
+    ASSERT_EQ(-3.0, L1.coeff(f,f));
+    ASSERT_EQ(0.0, L1.row(f).sum());
+    ASSERT_EQ(0.0, L1.col(f).sum());
+#else
+    //Soft assert if we have not edge_lenght_squared
+    ASSERT_NEAR(-3.0, L1.coeff(f,f), epsilon);
+    ASSERT_NEAR(0.0, L1.row(f).sum(), epsilon);
+    ASSERT_NEAR(0.0, L1.col(f).sum(), epsilon);
+#endif
+
+  }
+
+  //Same for huge cube.
+  igl::cotmatrix(V_huge,F,L1);
+  ASSERT_EQ(V.rows(), L1.rows());
+  ASSERT_EQ(V.rows(), L1.cols());
+  for(int f = 0;f<L1.rows();f++)
+  {
+    ASSERT_NEAR(-3.0, L1.coeff(f,f), epsilon);
+    ASSERT_NEAR(0.0, L1.row(f).sum(), epsilon);
+    ASSERT_NEAR(0.0, L1.col(f).sum(), epsilon);
+  }
+
+  //Same for tiny cube. we need to use a tolerance this time...
+  igl::cotmatrix(V_tiny,F,L1);
+  ASSERT_EQ(V.rows(), L1.rows());
+  ASSERT_EQ(V.rows(), L1.cols());
+  for(int f = 0;f<L1.rows();f++)
+  {
+    ASSERT_NEAR(-3.0, L1.coeff(f,f), epsilon);
+    ASSERT_NEAR(0.0, L1.row(f).sum(), epsilon);
+    ASSERT_NEAR(0.0, L1.col(f).sum(), epsilon);
+  }
+
+}//TEST(cotmatrix, cube)
+
+TEST(cotmatrix, tetrahedron)
+{
+  //The allowed error for this test
+  const double epsilon = 1e-15;
+
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  //This is a cube of dimensions 1.0x1.0x1.0
+  test_common::load_mesh("cube.obj", V, F);
+
+  //Prepare another mesh with triangles along side diagonals of the cube
+  //These triangles are form a regular tetrahedron of side sqrt(2)
+  Eigen::MatrixXi F_equi(4,3);
+  F_equi << 4,6,1,
+            6,4,3,
+            4,1,3,
+            1,6,3;
+
+  //Scale the cube to have huge sides
+  Eigen::MatrixXd V_huge = V * 1.0e8;
+
+  //Scale the cube to have tiny sides
+  Eigen::MatrixXd V_tiny = V * 1.0e-8;
+
+  //Check cotmatrix (Laplacian)
+  //The laplacian for the cube is quite singular.
+  //Each edge in a diagonal has two opposite angles of 90, with cotangent 0.0 each
+  //Each edge in a side has two opposite angle of 45, with (half)cotangen 0.5 each
+  //So the cotangent matrix always are (0+0) or (0.5+0.5)
+  Eigen::SparseMatrix<double> L1;
+
+  //Check the regular tetrahedron of side sqrt(2)
+  igl::cotmatrix(V,F_equi,L1);
+
+  ASSERT_EQ(V.rows(), L1.rows());
+  ASSERT_EQ(V.rows(), L1.cols());
+  for(int f = 0;f<L1.rows();f++)
+  {
+    //Check the diagonal. Only can value 0.0 for unused vertex or -3 / tan(60)
+    if (L1.coeff(f,f) < -0.1)
+        ASSERT_NEAR(-3 / tan(M_PI / 3.0), L1.coeff(f,f), epsilon);
+    else
+        ASSERT_NEAR(0.0, L1.coeff(f,f), epsilon);
+#ifdef IGL_EDGE_LENGTHS_SQUARED_H
+    //Hard assert if we have edge_lenght_squared
+    ASSERT_EQ(0.0, L1.row(f).sum());
+    ASSERT_EQ(0.0, L1.col(f).sum());
+#else
+    //Soft assert if we have not edge_lenght_squared
+    ASSERT_NEAR(0.0, L1.row(f).sum(), epsilon);
+    ASSERT_NEAR(0.0, L1.col(f).sum(), epsilon);
+#endif
+  }
+
+  //Check the huge regular tetrahedron
+  igl::cotmatrix(V_huge,F_equi,L1);
+
+  ASSERT_EQ(V.rows(), L1.rows());
+  ASSERT_EQ(V.rows(), L1.cols());
+  for(int f = 0;f<L1.rows();f++)
+  {
+    //Check the diagonal. Only can value 0.0 for unused vertex or -3 / tan(60)
+    if (L1.coeff(f,f) < -0.1)
+        ASSERT_NEAR(-3 / tan(M_PI / 3.0), L1.coeff(f,f), epsilon);
+    else
+        ASSERT_NEAR(0.0, L1.coeff(f,f), epsilon);
+    ASSERT_NEAR(0.0, L1.row(f).sum(), epsilon);
+    ASSERT_NEAR(0.0, L1.col(f).sum(), epsilon);
+  }
+
+  //Check the tiny regular tetrahedron
+  igl::cotmatrix(V_tiny,F_equi,L1);
+
+  ASSERT_EQ(V.rows(), L1.rows());
+  ASSERT_EQ(V.rows(), L1.cols());
+  for(int f = 0;f<L1.rows();f++)
+  {
+    //Check the diagonal. Only can value 0.0 for unused vertex or -3 / tan(60)
+    if (L1.coeff(f,f) < -0.1)
+        ASSERT_NEAR(-3 / tan(M_PI / 3.0), L1.coeff(f,f), epsilon);
+    else
+        ASSERT_NEAR(0.0, L1.coeff(f,f), epsilon);
+    ASSERT_NEAR(0.0, L1.row(f).sum(), epsilon);
+    ASSERT_NEAR(0.0, L1.col(f).sum(), epsilon);
+  }
+
+}//TEST(cotmatrix, tetrahedron)

tests/include/igl/cotmatrix_entries.cpp

@@ -0,0 +1,136 @@
+#include <test_common.h>
+#include <igl/cotmatrix_entries.h>
+
+TEST(cotmatrix_entries, simple)
+{
+  //The allowed error for this test
+  const double epsilon = 1e-15;
+
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  //This is a cube of dimensions 1.0x1.0x1.0
+  test_common::load_mesh("cube.obj", V, F);
+
+  //Prepare another mesh with triangles along side diagonals of the cube
+  //These triangles are form a regular tetrahedron of side sqrt(2)
+  Eigen::MatrixXi F_tet(4,3);
+  F_tet << 4,6,1,
+            6,4,3,
+            4,1,3,
+            1,6,3;
+
+  //1. Check cotmatrix_entries
+
+  Eigen::MatrixXd C1;
+  igl::cotmatrix_entries(V,F,C1);
+  ASSERT_EQ(F.rows(), C1.rows());
+  ASSERT_EQ(3, C1.cols());
+  //All angles in unit cube measure 45 or 90 degrees
+  //Their (half)cotangent must value 0.5 or 0.0
+  for(int f = 0;f<C1.rows();f++)
+  {
+#ifdef IGL_EDGE_LENGTHS_SQUARED_H
+      //Hard assert if we have edge_lenght_squared
+      for(int v = 0;v<3;v++)
+        if (C1(f,v) > 0.1)
+            ASSERT_EQ(0.5, C1(f,v));
+        else
+           ASSERT_EQ(0.0, C1(f,v));
+       //All cotangents sum 1.0 for those triangles
+       ASSERT_EQ(1.0, C1.row(f).sum());
+#else
+      //Soft assert if we have not edge_lenght_squared
+      for(int v = 0;v<3;v++)
+        if (C1(f,v) > 0.1)
+            ASSERT_NEAR(0.5, C1(f,v), epsilon);
+        else
+            ASSERT_NEAR(0.0, C1(f,v), epsilon);
+       //All cotangents sum 1.0 for those triangles
+       ASSERT_NEAR(1.0, C1.row(f).sum(), epsilon);
+#endif
+  }
+
+  //Check the regular tetrahedron
+  Eigen::MatrixXd C2;
+  igl::cotmatrix_entries(V,F_tet,C2);
+  ASSERT_EQ(F_tet.rows(), C2.rows());
+  ASSERT_EQ(3, C2.cols());
+  for(int f = 0;f<C2.rows();f++)
+  {
+    //Their (half)cotangent must value 0.5 / tan(M_PI / 3.0)
+    for(int v = 0;v<3;v++)
+       ASSERT_NEAR(0.5 / tan(M_PI / 3.0), C2(f,v), epsilon);
+  }
+
+  //Scale the cube to have huge sides
+  Eigen::MatrixXd V_huge = V * 1.0e8;
+  igl::cotmatrix_entries(V_huge,F,C1);
+  ASSERT_EQ(F.rows(), C1.rows());
+  ASSERT_EQ(3, C1.cols());
+  //All angles still measure 45 or 90 degrees
+  //Their (half)cotangent must value 0.5 or 0.0
+  for(int f = 0;f<C1.rows();f++)
+  {    
+#ifdef IGL_EDGE_LENGTHS_SQUARED_H
+      //Hard assert if we have edge_lenght_squared
+      for(int v = 0;v<3;v++)
+        if (C1(f,v) > 0.1)
+            ASSERT_EQ(0.5, C1(f,v));
+        else
+           ASSERT_EQ(0.0, C1(f,v));
+       //All cotangents sum 1.0 for those triangles
+       ASSERT_EQ(1.0, C1.row(f).sum());
+#else
+      //Soft assert if we have not edge_lenght_squared
+      for(int v = 0;v<3;v++)
+        if (C1(f,v) > 0.1)
+            ASSERT_NEAR(0.5, C1(f,v), epsilon);
+        else
+            ASSERT_NEAR(0.0, C1(f,v), epsilon);
+       //All cotangents sum 1.0 for those triangles
+       ASSERT_NEAR(1.0, C1.row(f).sum(), epsilon);
+#endif
+
+  }
+
+  //Check the huge regular tetrahedron
+  igl::cotmatrix_entries(V_huge,F_tet,C2);
+  ASSERT_EQ(F_tet.rows(), C2.rows());
+  ASSERT_EQ(3, C2.cols());
+  for(int f = 0;f<C2.rows();f++)
+  {
+    //Their (half)cotangent must value 0.5 / tan(M_PI / 3.0)
+    for(int v = 0;v<3;v++)
+       ASSERT_NEAR(0.5 / tan(M_PI / 3.0), C2(f,v), epsilon);
+  }
+
+  //Scale the cube to have tiny sides
+  Eigen::MatrixXd V_tiny = V * 1.0e-8;
+  igl::cotmatrix_entries(V_tiny,F,C1);
+  ASSERT_EQ(F.rows(), C1.rows());
+  ASSERT_EQ(3, C1.cols());
+  //All angles still measure 45 or 90 degrees
+  //Their (half)cotangent must value 0.5 or 0.0
+  for(int f = 0;f<C1.rows();f++)
+  {
+    for(int v = 0;v<3;v++)
+        if (C1(f,v) > 0.1)
+           ASSERT_NEAR(0.5, C1(f,v), epsilon);
+        else
+           ASSERT_NEAR(0.0, C1(f,v), epsilon);
+    //All cotangents sum 1.0 for those triangles
+    ASSERT_NEAR(1.0, C1.row(f).sum(), epsilon);
+  }
+
+  //Check the tiny regular tetrahedron
+  igl::cotmatrix_entries(V_tiny,F_tet,C2);
+  ASSERT_EQ(F_tet.rows(), C2.rows());
+  ASSERT_EQ(3, C2.cols());
+  for(int f = 0;f<C2.rows();f++)
+  {
+    //Their (half)cotangent must value 0.5 / tan(M_PI / 3.0)
+    for(int v = 0;v<3;v++)
+       ASSERT_NEAR(0.5 / tan(M_PI / 3.0), C2(f,v), epsilon);
+  }
+
+}//TEST(cotmatrix_entries, simple)

tests/include/igl/cut_to_disk.cpp

@@ -0,0 +1,126 @@
+#include <test_common.h>
+#include <igl/cut_mesh.h>
+#include <igl/cut_to_disk.h>
+#include <igl/edges.h>
+
+#include <array>
+#include <iostream>
+#include <vector>
+
+namespace cut_to_disk_helper {
+  template<typename DerivedV, typename DerivedF>
+  void assert_is_disk(
+      const Eigen::PlainObjectBase<DerivedV>& V,
+      const Eigen::PlainObjectBase<DerivedF>& F,
+      const std::vector<std::vector<int>>& cuts) {
+    using namespace igl;
+    std::set<std::array<int, 2>> cut_edges;
+    for (const auto& cut : cuts) {
+      const size_t cut_len = cut.size();
+      for (size_t i=0; i<cut_len-1; i++) {
+        std::array<int, 2> e{cut[i], cut[i+1]};
+        if (e[0] > e[1]) {
+          std::swap(e[0], e[1]);
+        }
+        cut_edges.insert(e);
+      }
+    }
+
+    const size_t num_faces = F.rows();
+    Eigen::MatrixXi cut_mask(num_faces, 3);
+    cut_mask.setZero();
+    for (size_t i=0; i<num_faces; i++) {
+      std::array<int, 2> e0{F(i, 0), F(i, 1)};
+      std::array<int, 2> e1{F(i, 1), F(i, 2)};
+      std::array<int, 2> e2{F(i, 2), F(i, 0)};
+      if (e0[0] > e0[1]) std::swap(e0[0], e0[1]);
+      if (e1[0] > e1[1]) std::swap(e1[0], e1[1]);
+      if (e2[0] > e2[1]) std::swap(e2[0], e2[1]);
+
+      if (cut_edges.find(e0) != cut_edges.end()) {
+        cut_mask(i, 0) = 1;
+      }
+      if (cut_edges.find(e1) != cut_edges.end()) {
+        cut_mask(i, 1) = 1;
+      }
+      if (cut_edges.find(e2) != cut_edges.end()) {
+        cut_mask(i, 2) = 1;
+      }
+    }
+
+    Eigen::MatrixXd V2;
+    Eigen::MatrixXi F2;
+    igl::cut_mesh(V, F, cut_mask, V2, F2);
+
+    Eigen::MatrixXi E2;
+    edges(F2, E2);
+    const auto euler = V2.rows() - E2.rows() + F2.rows();
+    ASSERT_TRUE(1 == euler || 2 == euler);
+  }
+}
+
+TEST(cut_to_disk, simple_tet) {
+  using namespace igl;
+  Eigen::MatrixXi F(4, 3);
+  F << 0, 2, 1,
+       0, 3, 2,
+       1, 2, 3,
+       0, 1, 3;
+  std::vector<std::vector<int>> cuts;
+  cut_to_disk(F, cuts);
+  ASSERT_EQ(0, cuts.size());
+}
+
+TEST(cut_to_disk, two_disconnected_tet) {
+  using namespace igl;
+  Eigen::MatrixXi F(8, 3);
+  F << 0, 2, 1,
+       0, 3, 2,
+       1, 2, 3,
+       0, 1, 3,
+       4, 6, 5,
+       4, 7, 6,
+       5, 6, 7,
+       4, 5, 7;
+  std::vector<std::vector<int>> cuts;
+  cut_to_disk(F, cuts);
+  ASSERT_EQ(0, cuts.size());
+}
+
+TEST(cut_to_disk, simple_square) {
+  using namespace igl;
+  Eigen::MatrixXi F(2, 3);
+  F << 0, 1, 2,
+       2, 1, 3;
+  std::vector<std::vector<int>> cuts;
+  cut_to_disk(F, cuts);
+  ASSERT_EQ(1, cuts.size());
+  ASSERT_EQ(5, cuts[0].size());
+  ASSERT_EQ(cuts[0][0], cuts[0][4]);
+}
+
+TEST(cut_to_disk, torus) {
+  using namespace igl;
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  test_common::load_mesh("TinyTorus.obj", V, F);
+
+  std::vector<std::vector<int>> cuts;
+  cut_to_disk(F, cuts);
+  ASSERT_EQ(2, cuts.size());
+
+  cut_to_disk_helper::assert_is_disk(V, F, cuts);
+}
+
+TEST(cut_to_disk, cube) {
+  using namespace igl;
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  test_common::load_mesh("cube.obj", V, F);
+
+  std::vector<std::vector<int>> cuts;
+  cut_to_disk(F, cuts);
+  ASSERT_EQ(0, cuts.size());
+
+  cut_to_disk_helper::assert_is_disk(V, F, cuts);
+}

tests/include/igl/decimate.cpp

@@ -0,0 +1,77 @@
+#include <test_common.h>
+#include <igl/decimate.h>
+#include <igl/sort.h>
+#include <igl/sortrows.h>
+#include <igl/normalize_row_lengths.h>
+#include <igl/slice.h>
+#include <igl/matlab_format.h>
+#include <iostream>
+
+class decimate : public ::testing::TestWithParam<std::string> {};
+
+TEST(decimate,hemisphere)
+{
+  // Load a hemisphere centered at the origin. For each original vertex compute
+  // its "perfect normal" (i.e., its position treated as unit vectors).
+  // Decimate the model and using the birth indices of the output vertices grab
+  // their original "perfect normals" and compare them to their current
+  // positions treated as unit vectors. If vertices have not moved much, then
+  // these should be similar (mostly this is checking if the birth indices are
+  // sane).
+  Eigen::MatrixXd V,U;
+  Eigen::MatrixXi F,G;
+  Eigen::VectorXi J,I;
+  // Load example mesh: GetParam() will be name of mesh file
+  test_common::load_mesh("hemisphere.obj", V, F);
+  // Perfect normals from positions
+  Eigen::MatrixXd NV = V.rowwise().normalized();
+  // Remove half of the faces
+  igl::decimate(V,F,F.rows()/2,U,G,J,I);
+  // Expect that all normals still point in same direction as original
+  Eigen::MatrixXd NU = U.rowwise().normalized();
+  Eigen::MatrixXd NVI;
+  igl::slice(NV,I,1,NVI);
+  ASSERT_EQ(NVI.rows(),NU.rows());
+  ASSERT_EQ(NVI.cols(),NU.cols());
+  // Dot product
+  Eigen::VectorXd D = (NU.array()*NVI.array()).rowwise().sum();
+  Eigen::VectorXd O = Eigen::VectorXd::Ones(D.rows());
+  // 0.2 chosen to succeed on 256 face hemisphere.obj reduced to 128 faces
+  test_common::assert_near(D,O,0.02);
+}
+
+TEST_P(decimate, closed)
+{
+  Eigen::MatrixXd V,U;
+  Eigen::MatrixXi F,G;
+  Eigen::VectorXi J;
+  // Load example mesh: GetParam() will be name of mesh file
+  test_common::load_mesh(GetParam(), V, F);
+  igl::decimate(V,F,0,U,G,J);
+  ASSERT_EQ(U.rows(),4);
+  ASSERT_EQ(G.rows(),4);
+  {
+    Eigen::MatrixXi I;
+    igl::sort(Eigen::MatrixXi(G),2,true,G,I);
+  }
+  {
+    Eigen::VectorXi I;
+    igl::sortrows(Eigen::MatrixXi(G),true,G,I);
+  }
+  // Tet with sorted faces
+  Eigen::MatrixXi T(4,3);
+  T<<
+    0,1,2,
+    0,1,3,
+    0,2,3,
+    1,2,3;
+  test_common::assert_eq(G,T);
+}
+
+INSTANTIATE_TEST_CASE_P
+(
+  closed_genus_0_meshes,
+  decimate,
+  ::testing::ValuesIn(test_common::closed_genus_0_meshes()),
+  test_common::string_test_name
+);

tests/include/igl/dirname.cpp

@@ -0,0 +1,47 @@
+#include <test_common.h>
+#include <igl/dirname.h>
+
+#include <string>
+#include <vector>
+#include <tuple>
+
+TEST(dirname, examples)
+{
+  const std::vector<
+    std::tuple<std::string,std::string> >
+    examples = 
+  {
+    std::make_tuple("/foo"                 ,"/"           ),
+    std::make_tuple("/foo/"                ,"/"           ),
+    std::make_tuple("/foo//"               ,"/"           ),
+    std::make_tuple("/foo/./"              ,"/foo"        ),
+    std::make_tuple("/foo/bar"             ,"/foo"        ),
+    std::make_tuple("/foo/bar."            ,"/foo"        ),
+    std::make_tuple("/foo/bar.txt"         ,"/foo"        ),
+    std::make_tuple("/foo/bar.txt.zip"     ,"/foo"        ),
+    std::make_tuple("/foo/bar.dir/"        ,"/foo"        ),
+    std::make_tuple("/foo/bar.dir/file"    ,"/foo/bar.dir"),
+    std::make_tuple("/foo/bar.dir/file.txt","/foo/bar.dir"),
+    std::make_tuple("."                   ,"."           ),
+    std::make_tuple("../foo"              ,".."           ),
+    std::make_tuple("./foo"              ,"."           ),
+    std::make_tuple("../foo/"             ,".."           ),
+    std::make_tuple("foo/"               ,"."           ),
+    std::make_tuple("foo//"               ,"."           ),
+    std::make_tuple("foo/./"              ,"foo"        ),
+    std::make_tuple("foo/bar"             ,"foo"        ),
+    std::make_tuple("foo/bar."            ,"foo"        ),
+    std::make_tuple("foo/bar.txt"         ,"foo"        ),
+    std::make_tuple("foo/bar.txt.zip"     ,"foo"        ),
+    std::make_tuple("foo/bar.dir/"        ,"foo"        ),
+    std::make_tuple("foo/bar.dir/file"    ,"foo/bar.dir"),
+    std::make_tuple("foo/bar.dir/file.txt","foo/bar.dir")
+  };
+  for(const auto & example : examples)
+  {
+    std::string d;
+    d = igl::dirname(std::get<0>(example));
+    ASSERT_EQ(std::get<1>(example),d);
+  }
+}
+

tests/include/igl/doublearea.cpp

@@ -0,0 +1,39 @@
+#include <test_common.h>
+#include <igl/doublearea.h>
+
+class doublearea : public ::testing::TestWithParam<std::string> {};
+
+TEST_P(doublearea, VF_vs_ABC )
+{
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  test_common::load_mesh(GetParam(), V, F);
+
+  // Check that computing double area with (V,F) is the same as computing
+  // double area with (V1,V2,V2)
+  Eigen::VectorXd A1,A2;
+  igl::doublearea(V,F,A1);
+  Eigen::MatrixXd A(F.rows(),3);
+  Eigen::MatrixXd B(F.rows(),3);
+  Eigen::MatrixXd C(F.rows(),3);
+  for(int f = 0;f<F.rows();f++)
+  {
+    A.row(f) = V.row(F(f,0));
+    B.row(f) = V.row(F(f,1));
+    C.row(f) = V.row(F(f,2));
+  }
+  igl::doublearea(A,B,C,A2);
+  ASSERT_EQ(A1.size(),A2.size());
+  for(int a = 0;a<A1.size();a++)
+  {
+    ASSERT_NEAR(A1(a),A2(a),1e-15);
+  }
+}
+
+INSTANTIATE_TEST_CASE_P
+(
+ all_meshes,
+ doublearea,
+ ::testing::ValuesIn(test_common::all_meshes()),
+ test_common::string_test_name
+);

tests/include/igl/edge_flaps.cpp

@@ -0,0 +1,46 @@
+#include <test_common.h>
+#include <igl/edge_flaps.h>
+
+class edge_flaps : public ::testing::TestWithParam<std::string> {};
+
+TEST_P(edge_flaps, verify)
+{
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  test_common::load_mesh(GetParam(), V, F);
+
+  Eigen::MatrixXi efE,efEF,efEI;
+  Eigen::VectorXi efEMAP;
+  igl::edge_flaps(F,efE,efEMAP,efEF,efEI);
+  ASSERT_EQ(efE.rows(),efEF.rows());
+  ASSERT_EQ(efE.cols(),2);
+  ASSERT_EQ(efE.cols(),efEF.cols());
+  // for each edge, make sure edge appears in face
+  for(int e = 0;e<efE.rows();e++)
+  {
+    for(int fe = 0;fe<2;fe++)
+    {
+      const int f = efEF(e,fe);
+      // index of corner
+      const int c = efEI(e,fe);
+      ASSERT_TRUE(f<F.rows());
+      // only check if not on boundary
+      if(f >= 0)
+      {
+        EXPECT_TRUE( 
+        // Either efE(e,[1 2]) = [i,j] appears after vertex c of face f
+          ((efE(e,0) == F(f,(c+1)%3)) && (efE(e,1) == F(f,(c+2)%3))) ||
+        // Or  efE(e,[2 1]) = [j,i] appears after vertex c of face f
+          ((efE(e,1) == F(f,(c+1)%3)) && (efE(e,0) == F(f,(c+2)%3))));
+      }
+    }
+  }
+}
+
+INSTANTIATE_TEST_CASE_P
+(
+ all_meshes,
+ edge_flaps,
+ ::testing::ValuesIn(test_common::all_meshes()),
+ test_common::string_test_name
+);

tests/include/igl/edge_lengths.cpp

@@ -0,0 +1,85 @@
+#include <test_common.h>
+#include <igl/edge_lengths.h>
+#include <iostream>
+
+TEST(edge_lengths, cube)
+{
+  //The allowed error for this test
+  const double epsilon = 1e-15;
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  //This is a cube of dimensions 1.0x1.0x1.0
+  test_common::load_mesh("cube.obj", V, F);
+  //Create scaled versions of the cube
+  double scale = 1.0;
+  double huge_scale = 1.0e8;
+  double tiny_scale = 1.0e-8;
+  Eigen::MatrixXd V_huge = V * huge_scale;
+  Eigen::MatrixXd V_tiny = V * tiny_scale;
+  //Prepare another mesh with triangles along side diagonals of the cube
+  //These triangles are form a regular tetrahedron of side sqrt(2)
+  Eigen::MatrixXi F_tet(4,3);
+  F_tet << 4,6,1,
+            6,4,3,
+            4,1,3,
+            1,6,3;
+
+  //2. Check edge_lengths
+  double side = 1.0;       //lenght of a side
+  double diag = sqrt(2.0); //lenght of a diagonal
+  Eigen::MatrixXd L;
+  igl::edge_lengths(V,F,L);
+  ASSERT_EQ(F.rows(), L.rows());
+  ASSERT_EQ(3, L.cols());
+  //All edges in unit cube measure 1.0 or sqrt(2) in diagonals
+  for(int f = 0;f<L.rows();f++)
+  {
+    //All edge_lengths_squared must be exactly "side" or "diag"
+    for(int e = 0;e<3;e++)
+        if (L(f,e) > 1.1*side)
+           ASSERT_EQ(diag, L(f,e));
+        else
+           ASSERT_EQ(side, L(f,e));
+    //All sides sum exactly side + side + diag
+    ASSERT_EQ(L.row(f).sum(), side + side + diag);
+  }
+
+  //Check the cube of huge sides
+  scale = 1.0e8;
+  side = scale;       //lenght of a side
+  diag = scale*sqrt(2.0); //lenght of a diagonal
+  igl::edge_lengths(V_huge,F,L);
+  ASSERT_EQ(F.rows(), L.rows());
+  ASSERT_EQ(3, L.cols());
+  for(int f = 0;f<L.rows();f++)
+  {
+    //All edge_lengths_squared must be exactly "side" or "diag"
+    for(int e = 0;e<3;e++)
+        if (L(f,e) > 1.1*side)
+           ASSERT_EQ(diag, L(f,e));
+        else
+           ASSERT_EQ(side, L(f,e));
+    //All sides sum exactly side + side + diag
+    ASSERT_NEAR(L.row(f).sum(), side + side + diag, epsilon);
+  }
+
+  //Check the cube of tiny sides
+  scale = 1.0e-8;
+  side = scale;       //lenght of a side
+  diag = scale*sqrt(2.0); //lenght of a diagonal
+  igl::edge_lengths(V_tiny,F,L);
+  ASSERT_EQ(F.rows(), L.rows());
+  ASSERT_EQ(3, L.cols());
+  for(int f = 0;f<L.rows();f++)
+  {
+    //All edge_lengths_squared must be exactly "side" or "diag"
+    for(int e = 0;e<3;e++)
+        if (L(f,e) > 1.1*side)
+           ASSERT_EQ(diag, L(f,e));
+        else
+           ASSERT_EQ(side, L(f,e));
+    //All sides sum exactly side + side + diag
+    ASSERT_EQ(L.row(f).sum(), side + side + diag);
+  }
+
+}

tests/include/igl/guess_extension.cpp

@@ -0,0 +1,30 @@
+
+#include <test_common.h>
+#include <igl/guess_extension.h>
+#include <igl/pathinfo.h>
+#include <cstdio>
+
+class guess_extension : public ::testing::TestWithParam<std::string> {};
+
+TEST_P(guess_extension, all_meshes)
+{
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  std::string path(test_common::data_path(GetParam()));
+  // Load example mesh: GetParam() will be name of mesh file
+  std::string d,b,e,f;
+  igl::pathinfo(path,d,b,e,f);
+  // Convert extension to lower case
+  std::transform(e.begin(), e.end(), e.begin(), ::tolower);
+  FILE * fp = fopen(path.c_str(),"r");
+  std::string guess = igl::guess_extension(fp);
+  ASSERT_EQ(guess,e);
+}
+
+INSTANTIATE_TEST_CASE_P
+(
+  all_meshes,
+  guess_extension,
+  ::testing::ValuesIn(test_common::all_meshes()),
+  test_common::string_test_name
+);

tests/include/igl/is_edge_manifold.cpp

@@ -0,0 +1,29 @@
+#include <test_common.h>
+#include <igl/is_edge_manifold.h>
+
+class is_edge_manifold : public ::testing::TestWithParam<std::string> {};
+
+TEST_P(is_edge_manifold, positive)
+{
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  test_common::load_mesh(GetParam(), V, F);
+  ASSERT_TRUE( igl::is_edge_manifold(F) );
+}
+
+TEST(is_edge_manifold, negative)
+{
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  // Known non-manifold mesh
+  test_common::load_mesh("truck.obj", V, F);
+  ASSERT_FALSE( igl::is_edge_manifold(F) );
+}
+
+INSTANTIATE_TEST_CASE_P
+(
+ manifold_meshes,
+ is_edge_manifold,
+ ::testing::ValuesIn(test_common::manifold_meshes()),
+ test_common::string_test_name
+);

tests/include/igl/is_symmetric.cpp

@@ -0,0 +1,35 @@
+
+#include <test_common.h>
+#include <igl/is_symmetric.h>
+
+TEST(is_symmetric, sparse)
+{
+  {
+    Eigen::MatrixXd M(3,3);
+    M<<1,2,3,4,5,6,7,8,9;
+    Eigen::SparseMatrix<double> S = M.sparseView();
+    ASSERT_FALSE(igl::is_symmetric(S));
+  }
+  {
+    Eigen::MatrixXd M(3,3);
+    M<<1,2,3,2,4,5,3,5,6;
+    Eigen::SparseMatrix<double> S = M.sparseView();
+    ASSERT_FALSE(igl::is_symmetric(S));
+  }
+}
+
+TEST(is_symmetric, dense)
+{
+  {
+    Eigen::MatrixXd M(3,3);
+    M<<1,2,3,4,5,6,7,8,9;
+    ASSERT_FALSE(igl::is_symmetric(M));
+  }
+  {
+    Eigen::MatrixXd M(3,3);
+    M<<1,2,3,
+       2,4,5,
+       3,5,6;
+    ASSERT_FALSE(igl::is_symmetric(M));
+  }
+}

tests/include/igl/ismember.cpp

@@ -0,0 +1,44 @@
+#include <test_common.h>
+#include <igl/ismember.h>
+#include <igl/matlab_format.h>
+
+TEST(ismember, simple)
+{
+  Eigen::MatrixXi A(3,4);
+  A<<11,12,13,14,21,22,23,24,31,32,33,34;
+  Eigen::MatrixXi B(2,3);
+  B<<11,13,11,21,22,34;
+
+  Eigen::Matrix<bool,Eigen::Dynamic,Eigen::Dynamic> IA;
+  Eigen::MatrixXi LOCB;
+  igl::ismember(A,B,IA,LOCB);
+  Eigen::Map<Eigen::VectorXi> vB = 
+    Eigen::Map<Eigen::VectorXi>(B.data(),B.rows()*B.cols());
+  for(int i = 0;i<A.rows();i++)
+  {
+    for(int j = 0;j<A.cols();j++)
+    {
+      // try to find in b
+      int bi = 0;
+      for(;bi<vB.size();bi++)
+      {
+        if(A(i,j) == vB(bi))
+        {
+          break;
+        }
+      }
+      if(IA(i,j))
+      {
+        ASSERT_GE(LOCB(i,j),0);
+        ASSERT_LT(LOCB(i,j),B.size());
+        ASSERT_EQ(vB(LOCB(i,j)),A(i,j));
+        ASSERT_EQ(LOCB(i,j),bi);
+      }else
+      {
+        ASSERT_EQ(LOCB(i,j),-1);
+        ASSERT_EQ(bi,B.size());
+      }
+    }
+  }
+}
+

tests/include/igl/list_to_matrix.cpp

@@ -0,0 +1,54 @@
+#include <test_common.h>
+#include <igl/list_to_matrix.h>
+#include <igl/STR.h>
+#include <tuple>
+
+namespace list_to_matrix
+{
+  typedef std::tuple<int/*n*/,int/*m*/> NM;
+  inline std::string NM_test_name(
+    const ::testing::TestParamInfo<NM>& info)
+  {
+    return STR(
+      std::get<0>(info.param)<<"x"<<
+      std::get<1>(info.param)<<"_");
+  };
+}
+class ListToMatrixTest : public ::testing::TestWithParam<list_to_matrix::NM> {};
+
+TEST_P(ListToMatrixTest,matrix)
+{
+  const int n = std::get<0>(GetParam());
+  const int m = std::get<1>(GetParam());
+  std::vector<std::vector<double> > vX(n,std::vector<double>(m));
+  for(int i = 0;i<n;i++)
+  {
+    for(int j = 0;j<m;j++)
+    {
+      vX[i][j] = i+j*n;
+    }
+  }
+  Eigen::MatrixXd mX;
+  igl::list_to_matrix(vX,mX);
+  for(int i = 0;i<n;i++)
+  {
+    for(int j = 0;j<m;j++)
+    {
+      ASSERT_EQ(vX[i][j],mX(i,j));
+    }
+  }
+}
+
+INSTANTIATE_TEST_CASE_P
+(
+  suite,
+  ListToMatrixTest,
+  ::testing::ValuesIn<std::vector<list_to_matrix::NM> >(
+    std::vector<list_to_matrix::NM>{
+    list_to_matrix::NM{100,4},
+    list_to_matrix::NM{4,100},
+    list_to_matrix::NM{100,1},
+    list_to_matrix::NM{1,100},
+    }),
+  list_to_matrix::NM_test_name
+);

tests/include/igl/main.cpp

@@ -0,0 +1,6 @@
+#include <gtest/gtest.h>
+
+int main(int argc, char **argv) {
+      ::testing::InitGoogleTest(&argc, argv);
+        return RUN_ALL_TESTS();
+}

tests/include/igl/mosek/CMakeLists.txt

@@ -0,0 +1,6 @@
+file(GLOB TEST_SRC_FILES *.cpp main.cpp)
+file(GLOB TEST_INC_FILES *.h *.inl)
+
+add_executable(igl_mosek_tests ${TEST_SRC_FILES} ${TEST_INC_FILES})
+target_link_libraries(igl_mosek_tests igl::core igl::mosek gtest_main)
+add_test(NAME run_igl_mosek_tests COMMAND igl_mosek_tests)

tests/include/igl/mosek/bbw.cpp

@@ -0,0 +1,26 @@
+#include <test_common.h>
+#include <igl/boundary_conditions.h>
+#include <igl/readMESH.h>
+#include <igl/writeDMAT.h>
+#include <igl/readTGF.h>
+#include <igl/mosek/bbw.h>
+
+TEST(mosek_bbw, decimated_knight)
+{
+  Eigen::MatrixXd V,C;
+  Eigen::MatrixXi T,F,E;
+  igl::readMESH(test_common::data_path("decimated-knight.mesh"),V,T,F);
+  igl::readTGF(test_common::data_path("decimated-knight.tgf"),C,E);
+  Eigen::MatrixXd W_groundtruth, Was, Wmo;
+  igl::readDMAT(
+    test_common::data_path("decimated-knight-matlab-active-set.dmat"),W_groundtruth);
+  Eigen::VectorXi b;
+  Eigen::MatrixXd bc;
+  igl::boundary_conditions(V,T,C,Eigen::VectorXi(),E,Eigen::MatrixXi(),b,bc);
+  igl::BBWData params;
+  igl::mosek::MosekData mosek_params;
+  igl::mosek::bbw(V,T,b,bc,params,mosek_params,Wmo);
+  igl::writeDMAT("decimated-knight-mo.dmat",Wmo);
+  // Mosek is less accurate
+  ASSERT_LT( (Wmo-W_groundtruth).array().abs().maxCoeff() ,1e-3);
+}

tests/include/igl/mosek/main.cpp

@@ -0,0 +1,6 @@
+#include <gtest/gtest.h>
+
+int main(int argc, char **argv) {
+      ::testing::InitGoogleTest(&argc, argv);
+        return RUN_ALL_TESTS();
+}

tests/include/igl/pathinfo.cpp

@@ -0,0 +1,35 @@
+#include <test_common.h>
+#include <igl/pathinfo.h>
+
+#include <string>
+#include <vector>
+#include <tuple>
+
+TEST(pathinfo, examples)
+{
+  const std::vector<
+    std::tuple<std::string,std::string,std::string,std::string,std::string> >
+    examples = 
+  {
+    std::make_tuple("/foo"                 ,"/"           ,"foo"        ,""   ,"foo"),
+    std::make_tuple("/foo/"                ,"/"           ,"foo"        ,""   ,"foo"),
+    std::make_tuple("/foo//"               ,"/"           ,"foo"        ,""   ,"foo"),
+    std::make_tuple("/foo/./"              ,"/foo"        ,"."          ,""   ,""),
+    std::make_tuple("/foo/bar"             ,"/foo"        ,"bar"        ,""   ,"bar"),
+    std::make_tuple("/foo/bar."            ,"/foo"        ,"bar."       ,""   ,"bar"),
+    std::make_tuple("/foo/bar.txt"         ,"/foo"        ,"bar.txt"    ,"txt","bar"),
+    std::make_tuple("/foo/bar.txt.zip"     ,"/foo"        ,"bar.txt.zip","zip","bar.txt"),
+    std::make_tuple("/foo/bar.dir/"        ,"/foo"        ,"bar.dir"    ,"dir","bar"),
+    std::make_tuple("/foo/bar.dir/file"    ,"/foo/bar.dir","file"       ,""   ,"file"),
+    std::make_tuple("/foo/bar.dir/file.txt","/foo/bar.dir","file.txt"   ,"txt","file")
+  };
+  for(const auto & example : examples)
+  {
+    std::string d,b,e,f;
+    igl::pathinfo(std::get<0>(example),d,b,e,f);
+    ASSERT_EQ(std::get<1>(example),d);
+    ASSERT_EQ(std::get<2>(example),b);
+    ASSERT_EQ(std::get<3>(example),e);
+    ASSERT_EQ(std::get<4>(example),f);
+  }
+}

tests/include/igl/per_face_normals.cpp

@@ -0,0 +1,38 @@
+
+#include <test_common.h>
+#include <igl/per_face_normals.h>
+#include <Eigen/Geometry>
+
+class per_face_normals : public ::testing::TestWithParam<std::string> {};
+
+TEST_P(per_face_normals, dot)
+{
+  Eigen::MatrixXd V,N;
+  Eigen::MatrixXi F;
+  // Load example mesh: GetParam() will be name of mesh file
+  test_common::load_mesh(GetParam(), V, F);
+  igl::per_face_normals(V,F,N);
+  ASSERT_EQ(F.rows(),N.rows());
+  for(int f = 0;f<N.rows();f++)
+  {
+    for(int c = 0;c<3;c++)
+    {
+      // Every half-edge dot the normal should be 0
+      ASSERT_LT(
+        std::abs((V.row(F(f,c))-V.row(F(f,(c+1)%3))).dot(N.row(f))),
+        1e-12);
+    }
+  }
+  // ASSERT_EQ(a,b);
+  // ASSERT_TRUE(a==b);
+  // ASSERT_NEAR(a,b,1e-15)
+  // ASSERT_LT(a,1e-12);
+}
+
+INSTANTIATE_TEST_CASE_P
+(
+  all_meshes,
+  per_face_normals,
+  ::testing::ValuesIn(test_common::all_meshes()),
+  test_common::string_test_name
+);

tests/include/igl/qslim.cpp

@@ -0,0 +1,44 @@
+#include <test_common.h>
+#include <igl/qslim.h>
+#include <igl/cylinder.h>
+#include <igl/upsample.h>
+#include <igl/point_mesh_squared_distance.h>
+//#include <igl/hausdorff.h>
+#include <igl/writePLY.h>
+
+TEST(qslim,cylinder)
+{
+  using namespace igl;
+  const int axis_devisions = 5;
+  const int height_devisions = 2+10;
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  cylinder(axis_devisions,height_devisions,V,F);
+  Eigen::MatrixXd U;
+  Eigen::MatrixXi G;
+  Eigen::VectorXi I,J;
+  qslim(V,F,2*axis_devisions,U,G,I,J);
+  ASSERT_EQ(axis_devisions*2,U.rows());
+  double l,u;
+  igl::writePLY("qslim-cylinder-vf.ply",V,F);
+  igl::writePLY("qslim-cylinder-ug.ply",U,G);
+  const auto & hausdorff_lower_bound = [](
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & F,
+    Eigen::MatrixXd & U,
+    Eigen::MatrixXi & G)->double
+  {
+    Eigen::MatrixXd Vk;
+    Eigen::MatrixXi Fk;
+    igl::upsample(V,F,Vk,Fk,5);
+    Eigen::MatrixXd C;
+    Eigen::VectorXi I;
+    Eigen::VectorXd D;
+    igl::point_mesh_squared_distance(Vk,U,G,D,I,C);
+    return D.array().sqrt().maxCoeff();
+  };
+  //igl::hausdorff(V,F,U,G,1e-14,l,u);
+  ASSERT_NEAR(hausdorff_lower_bound(V,F,U,G),0,2e-10);
+  //igl::hausdorff(U,G,V,F,1e-14,l,u);
+  ASSERT_NEAR(hausdorff_lower_bound(U,G,V,F),0,2e-10);
+}

tests/include/igl/readDMAT.cpp

@@ -0,0 +1,19 @@
+#include <test_common.h>
+
+TEST(readDMAT, Comp) {
+    Eigen::MatrixXd N1, N2;
+    test_common::load_matrix("duplicated_faces_N1.dmat", N1);
+    test_common::load_matrix("duplicated_faces_N2.dmat", N2);
+
+    ASSERT_EQ(N1.rows(), N2.rows());
+    ASSERT_EQ(N1.cols(), N2.cols());
+    ASSERT_FALSE(((N1-N2).array() != 0.0).all());
+
+    const size_t rows = N1.rows();
+    const size_t cols = N1.cols();
+    for (size_t i=0; i<rows; i++) {
+        for (size_t j=0; j<cols; j++) {
+            ASSERT_FLOAT_EQ(N1(i,j), N2(i,j));
+        }
+    }
+}

tests/include/igl/readOBJ.cpp

@@ -0,0 +1,11 @@
+#include <test_common.h>
+
+TEST(readOBJ, simple) {
+    Eigen::MatrixXd V;
+    Eigen::MatrixXi F;
+    // wait... so this is actually testing test_common::load_mesh not readOBJ
+    // directly...
+    test_common::load_mesh("cube.obj", V, F);
+    ASSERT_EQ(8, V.rows());
+    ASSERT_EQ(12, F.rows());
+}

tests/include/igl/readOFF.cpp

@@ -0,0 +1,13 @@
+#include <test_common.h>
+
+TEST(readOFF, simple) {
+    Eigen::MatrixXd V;
+    Eigen::MatrixXi F;
+    // wait... so this is actually testing test_common::load_mesh not readOFF
+    // directly...
+    test_common::load_mesh("cube.off", V, F);
+    ASSERT_EQ(8, V.rows());   //has 8 vertex
+    ASSERT_EQ(3, V.cols());   //3D coordinates
+    ASSERT_EQ(12, F.rows());  //has 6*2=12 facets
+    ASSERT_EQ(3, F.cols());   //facets are triangles
+}

tests/include/igl/seam_edges.cpp

@@ -0,0 +1,23 @@
+#include <test_common.h>
+#include <igl/seam_edges.h>
+#include <igl/readOBJ.h>
+
+TEST(seam_edges, tet)
+{
+  Eigen::MatrixXd V,TC,CN;
+  Eigen::MatrixXi F,FTC,FN;
+  // Load example mesh: GetParam() will be name of mesh file
+  igl::readOBJ(test_common::data_path("tet.obj"), V, TC,CN,F,FTC,FN);
+  Eigen::MatrixXi seams,boundaries,foldovers;
+  igl::seam_edges(V,TC,F,FTC,seams,boundaries,foldovers);
+
+  Eigen::MatrixXi seams_gt(3,4);
+  seams_gt<<
+    0,0,1,2,
+    3,0,0,2,
+    1,0,3,2;
+  test_common::assert_eq(seams,seams_gt);
+  ASSERT_EQ(boundaries.size(),0);
+  ASSERT_EQ(foldovers.size(),0);
+}
+

tests/include/igl/setdiff.cpp

@@ -0,0 +1,46 @@
+#include <test_common.h>
+#include <igl/setdiff.h>
+
+TEST(setdiff,matrix)
+{
+  // Base cases
+  {
+    const Eigen::VectorXi A = (Eigen::VectorXi(4)<<1,2,1,3).finished();
+    const Eigen::VectorXi B(0,1);
+    Eigen::VectorXi C,IA;
+    const Eigen::VectorXi cC = (Eigen::VectorXi(3)<<1,2,3).finished();
+    const Eigen::VectorXi cIA = (Eigen::VectorXi(3)<<0,1,3).finished();
+    igl::setdiff(A,B,C,IA);
+    test_common::assert_eq(C,cC);
+    test_common::assert_eq(IA,cIA);
+  }
+  {
+    const Eigen::VectorXi A(0,1);
+    const Eigen::VectorXi B = (Eigen::VectorXi(4)<<1,2,1,3).finished();
+    Eigen::VectorXi C,IA;
+    const Eigen::VectorXi cC(0,1);
+    const Eigen::VectorXi cIA(0,1);
+    igl::setdiff(A,B,C,IA);
+    test_common::assert_eq(C,cC);
+    test_common::assert_eq(IA,cIA);
+  }
+
+  {
+    // Monkey test
+    Eigen::VectorXi A(12);
+    A = (Eigen::VectorXd::Random(A.size(),1).array().abs()*9).cast<int>();
+    Eigen::VectorXi B(12);
+    B = (Eigen::VectorXd::Random(B.size(),1).array().abs()*9).cast<int>();
+    Eigen::VectorXi C,IA;
+    igl::setdiff(A,B,C,IA);
+    for(int i = 0;i<C.size();i++)
+    {
+      ASSERT_EQ(C(i),A(IA(i)));
+      for(int j = 0;j<B.size();j++)
+      {
+        ASSERT_NE(B(j),C(i));
+      }
+    }
+  }
+}
+

tests/include/igl/slice.cpp

@@ -0,0 +1,109 @@
+#include <test_common.h>
+#include <igl/slice.h>
+#include <igl/LinSpaced.h>
+
+TEST(slice, dense_identity)
+{
+  // https://en.wikipedia.org/wiki/Monkey_testing
+  Eigen::MatrixXd A = Eigen::MatrixXd::Random(10,9);
+  Eigen::VectorXi I = igl::LinSpaced<Eigen::VectorXi >(A.rows(),0,A.rows()-1);
+  Eigen::VectorXi J = igl::LinSpaced<Eigen::VectorXi >(A.cols(),0,A.cols()-1);
+  {
+    Eigen::MatrixXd B;
+    igl::slice(A,I,J,B);
+    test_common::assert_eq(A,B);
+  }
+  {
+    Eigen::MatrixXd B;
+    igl::slice(A,I,1,B);
+    test_common::assert_eq(A,B);
+  }
+  {
+    Eigen::MatrixXd B;
+    igl::slice(A,J,2,B);
+    test_common::assert_eq(A,B);
+  }
+}
+
+TEST(slice, sparse_identity)
+{
+  Eigen::SparseMatrix<double> A = Eigen::MatrixXd::Random(10,9).sparseView();
+  Eigen::VectorXi I = igl::LinSpaced<Eigen::VectorXi >(A.rows(),0,A.rows()-1);
+  Eigen::VectorXi J = igl::LinSpaced<Eigen::VectorXi >(A.cols(),0,A.cols()-1);
+  {
+    Eigen::SparseMatrix<double> B;
+    igl::slice(A,I,J,B);
+    test_common::assert_eq(A,B);
+  }
+  {
+    Eigen::SparseMatrix<double> B;
+    igl::slice(A,I,1,B);
+    test_common::assert_eq(A,B);
+  }
+  {
+    Eigen::SparseMatrix<double> B;
+    igl::slice(A,J,2,B);
+    test_common::assert_eq(A,B);
+  }
+}
+
+TEST(slice,density_reverse)
+{
+  {
+    Eigen::MatrixXd A = Eigen::MatrixXd::Random(10,9);
+    Eigen::VectorXi I = igl::LinSpaced<Eigen::VectorXi >(A.rows(),A.rows()-1,0);
+    Eigen::VectorXi J = igl::LinSpaced<Eigen::VectorXi >(A.cols(),0,A.cols()-1);
+    Eigen::MatrixXd B;
+    igl::slice(A,I,J,B);
+    // reverse rows (i.e., reverse each column vector)
+    Eigen::MatrixXd C = A.colwise().reverse().eval();
+    test_common::assert_eq(B,C);
+  }
+  {
+    Eigen::MatrixXd A = Eigen::MatrixXd::Random(10,9);
+    Eigen::VectorXi I = igl::LinSpaced<Eigen::VectorXi >(A.rows(),0,A.rows()-1);
+    Eigen::VectorXi J = igl::LinSpaced<Eigen::VectorXi >(A.cols(),A.cols()-1,0);
+    Eigen::MatrixXd B;
+    igl::slice(A,I,J,B);
+    // reverse cols (i.e., reverse each row vector)
+    Eigen::MatrixXd C = A.rowwise().reverse().eval();
+    test_common::assert_eq(B,C);
+  }
+}
+
+TEST(slice,random)
+{
+  // Test whether unsorted indices are handled correctly by Randomly grow and
+  // shrink a matrix by slicing out rows and columns: note that growing will
+  // test whether repeated indices are correctly handled
+  std::vector<std::pair<int,int> > sizes = {{30,27},{3,4}};
+  for(const auto & size : sizes)
+  {
+    Eigen::MatrixXd A(10,9);
+    for(int i = 0;i<A.rows();i++)
+    {
+      for(int j = 0;j<A.cols();j++)
+      {
+        A(i,j) = A.rows()*j + i;
+      }
+    }
+    Eigen::VectorXi I = 
+      ((Eigen::VectorXd::Random(size.first,1).array()*0.5+0.5)*A.rows()
+       ).cast<int>();
+    Eigen::VectorXi J = 
+      ((Eigen::VectorXd::Random(size.second,1).array()*0.5+0.5)*A.cols()
+       ).cast<int>();
+    Eigen::MatrixXd B;
+    igl::slice(A,I,J,B);
+    Eigen::MatrixXd C(I.size(),J.size());
+    for(int i = 0;i<I.size();i++)
+    {
+      for(int j = 0;j<J.size();j++)
+      {
+        C(i,j) = A.rows()*J(j) + I(i);
+      }
+    }
+    test_common::assert_eq(B,C);
+  }
+}
+

tests/include/igl/slice_into.cpp

@@ -0,0 +1,72 @@
+#include <test_common.h>
+#include <igl/slice_into.h>
+#include <igl/LinSpaced.h>
+
+TEST(slice_into, dense_identity)
+{
+  Eigen::MatrixXd A = Eigen::MatrixXd::Random(10,9);
+  Eigen::VectorXi I = igl::LinSpaced<Eigen::VectorXi >(A.rows(),0,A.rows()-1);
+  Eigen::VectorXi J = igl::LinSpaced<Eigen::VectorXi >(A.cols(),0,A.cols()-1);
+  {
+    Eigen::MatrixXd B(I.maxCoeff()+1,J.maxCoeff()+1);
+    igl::slice_into(A,I,J,B);
+    test_common::assert_eq(A,B);
+  }
+  {
+    Eigen::MatrixXd B(I.maxCoeff()+1,A.cols());
+    igl::slice_into(A,I,1,B);
+    test_common::assert_eq(A,B);
+  }
+  {
+    Eigen::MatrixXd B(A.rows(),J.maxCoeff()+1);
+    igl::slice_into(A,J,2,B);
+    test_common::assert_eq(A,B);
+  }
+}
+
+TEST(slice_into,density_reverse)
+{
+  {
+    Eigen::MatrixXd A = Eigen::MatrixXd::Random(10,9);
+    Eigen::VectorXi I = igl::LinSpaced<Eigen::VectorXi >(A.rows(),A.rows()-1,0);
+    Eigen::VectorXi J = igl::LinSpaced<Eigen::VectorXi >(A.cols(),0,A.cols()-1);
+    Eigen::MatrixXd B(I.maxCoeff()+1,J.maxCoeff()+1);
+    igl::slice_into(A,I,J,B);
+    // reverse rows (i.e., reverse each column vector)
+    Eigen::MatrixXd C = A.colwise().reverse().eval();
+    test_common::assert_eq(B,C);
+  }
+  {
+    Eigen::MatrixXd A = Eigen::MatrixXd::Random(10,9);
+    Eigen::VectorXi I = igl::LinSpaced<Eigen::VectorXi >(A.rows(),0,A.rows()-1);
+    Eigen::VectorXi J = igl::LinSpaced<Eigen::VectorXi >(A.cols(),A.cols()-1,0);
+    Eigen::MatrixXd B(I.maxCoeff()+1,J.maxCoeff()+1);
+    igl::slice_into(A,I,J,B);
+    // reverse cols (i.e., reverse each row vector)
+    Eigen::MatrixXd C = A.rowwise().reverse().eval();
+    test_common::assert_eq(B,C);
+  }
+}
+
+
+TEST(slice_into, sparse_identity)
+{
+  Eigen::SparseMatrix<double> A = Eigen::MatrixXd::Random(10,9).sparseView();
+  Eigen::VectorXi I = igl::LinSpaced<Eigen::VectorXi >(A.rows(),0,A.rows()-1);
+  Eigen::VectorXi J = igl::LinSpaced<Eigen::VectorXi >(A.cols(),0,A.cols()-1);
+  {
+    Eigen::SparseMatrix<double> B(I.maxCoeff()+1,J.maxCoeff()+1);
+    igl::slice_into(A,I,J,B);
+    test_common::assert_eq(A,B);
+  }
+  {
+    Eigen::SparseMatrix<double> B(I.maxCoeff()+1,A.cols());
+    igl::slice_into(A,I,1,B);
+    test_common::assert_eq(A,B);
+  }
+  {
+    Eigen::SparseMatrix<double> B(A.rows(),J.maxCoeff()+1);
+    igl::slice_into(A,J,2,B);
+    test_common::assert_eq(A,B);
+  }
+}

tests/include/igl/sort.cpp

@@ -0,0 +1,91 @@
+#include <test_common.h>
+#include <igl/sort.h>
+#include <igl/STR.h>
+#include <tuple>
+
+namespace sort
+{
+  typedef std::tuple<int/*n*/,int/*m*/,int/*dim*/,bool/*ascending*/>
+    NMDimAscending;
+  inline std::string NMDimAscending_test_name(
+    const ::testing::TestParamInfo<NMDimAscending>& info)
+  {
+    return STR(
+      std::get<0>(info.param)<<"x"<<
+      std::get<1>(info.param)<<"_"<<
+      "dim_"<<std::get<2>(info.param)<<"_"<<
+      "ascending_"<<(std::get<3>(info.param)?"true":"false"));
+  };
+}
+class SortTest : public ::testing::TestWithParam<sort::NMDimAscending> {};
+
+TEST_P(SortTest,random)
+{
+  const int n = std::get<0>(GetParam());
+  const int m = std::get<1>(GetParam());
+  const int dim = std::get<2>(GetParam());
+  const bool ascending = std::get<3>(GetParam());
+  Eigen::MatrixXd X = Eigen::MatrixXd::Random(n,m);
+  // sort ascending
+  Eigen::MatrixXd Y;
+  Eigen::MatrixXi IX;
+  igl::sort(X,dim,ascending,Y,IX);
+  ASSERT_EQ(X.rows(),Y.rows());
+  ASSERT_EQ(X.cols(),Y.cols());
+  ASSERT_EQ(X.rows(),IX.rows());
+  ASSERT_EQ(X.cols(),IX.cols());
+  for(int i = 0;i<n;i++)
+  {
+    for(int j = 0;j<m;j++)
+    {
+      ASSERT_EQ(Y(i,j),X(dim==1?IX(i,j):i,dim==2?IX(i,j):j));
+    }
+  }
+  for(int i = (dim==1?1:0);i<n;i++) 
+  {
+    for(int j = (dim==2?1:0);j<m;j++)
+    {
+      if(ascending)
+      {
+        ASSERT_LE(Y(i-(dim==1?1:0),j-(dim==2?1:0)),Y(i,j));
+      }else
+      {
+        ASSERT_GE(Y(i-(dim==1?1:0),j-(dim==2?1:0)),Y(i,j));
+      }
+    }
+  }
+}
+
+INSTANTIATE_TEST_CASE_P
+(
+  suite,
+  SortTest,
+  ::testing::ValuesIn<std::vector<sort::NMDimAscending> >(
+    std::vector<sort::NMDimAscending> {
+    sort::NMDimAscending{100,3,1,true},
+    sort::NMDimAscending{100,3,2,true},
+    sort::NMDimAscending{100,3,1,false},
+    sort::NMDimAscending{100,3,2,false},
+    sort::NMDimAscending{3,100,1,true},
+    sort::NMDimAscending{3,100,2,true},
+    sort::NMDimAscending{3,100,1,false},
+    sort::NMDimAscending{3,100,2,false},
+    sort::NMDimAscending{100,2,1,true},
+    sort::NMDimAscending{100,2,2,true},
+    sort::NMDimAscending{100,2,1,false},
+    sort::NMDimAscending{100,2,2,false},
+    sort::NMDimAscending{2,100,1,true},
+    sort::NMDimAscending{2,100,2,true},
+    sort::NMDimAscending{2,100,1,false},
+    sort::NMDimAscending{2,100,2,false},
+    sort::NMDimAscending{100,4,1,true},
+    sort::NMDimAscending{100,4,2,true},
+    sort::NMDimAscending{100,4,1,false},
+    sort::NMDimAscending{100,4,2,false},
+    sort::NMDimAscending{4,100,1,true},
+    sort::NMDimAscending{4,100,2,true},
+    sort::NMDimAscending{4,100,1,false},
+    sort::NMDimAscending{4,100,2,false},
+    }),
+  sort::NMDimAscending_test_name
+);

tests/include/igl/squared_edge_lengths.cpp

@@ -0,0 +1,121 @@
+#include <test_common.h>
+#include <igl/squared_edge_lengths.h>
+#include <iostream>
+
+
+TEST(squared_edge_lengths, cube)
+{
+  //The allowed error for this test
+  const double epsilon = 1e-15;
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  //This is a cube of dimensions 1.0x1.0x1.0
+  test_common::load_mesh("cube.obj", V, F);
+  //Create scaled versions of the cube
+  double scale = 1.0;
+  double huge_scale = 1.0e8;
+  double tiny_scale = 1.0e-8;
+  Eigen::MatrixXd V_huge = V * huge_scale;
+  Eigen::MatrixXd V_tiny = V * tiny_scale;
+  //Prepare another mesh with triangles along side diagonals of the cube
+  //These triangles are form a regular tetrahedron of side sqrt(2)
+  Eigen::MatrixXi F_tet(4,3);
+  F_tet << 4,6,1,
+            6,4,3,
+            4,1,3,
+            1,6,3;
+
+  //1. Check edge_lengths_squared function
+  double side_sq = 1.0; //squared lenght of a side
+  double diag_sq = 2.0; //squared lenght of a diagonal
+  Eigen::MatrixXd L_sq;
+  igl::squared_edge_lengths(V,F,L_sq);
+  ASSERT_EQ(F.rows(), L_sq.rows());
+  ASSERT_EQ(3, L_sq.cols());
+  //All edges in unit cube measure 1.0 or sqrt(2) in diagonals
+  for(int f = 0;f<L_sq.rows();f++)
+  {
+    //All edge_lengths_squared must be exactly "side_sq" or "diag_sq"
+    for(int e = 0;e<3;e++)
+        if (L_sq(f,e) > 1.1)
+           ASSERT_EQ(diag_sq, L_sq(f,e));
+        else
+           ASSERT_EQ(side_sq, L_sq(f,e));
+    //All sides sum exactly side_sq + side_sq + diag_sq
+    ASSERT_EQ(L_sq.row(f).sum(), side_sq + side_sq + diag_sq);
+  }
+
+  //Check the regular tetrahedron
+  igl::squared_edge_lengths(V,F_tet,L_sq);
+  ASSERT_EQ(F_tet.rows(), L_sq.rows());
+  ASSERT_EQ(3, L_sq.cols());
+  //All edges measure sqrt(2)
+  for(int f = 0;f<L_sq.rows();f++)
+  {
+      //All edge_lengths_squared must be exactly "diag_sq"
+    for(int e = 0;e<3;e++)
+       ASSERT_EQ(2.0, L_sq(f,e));
+  }
+
+
+  //Scale the cube to have huge sides
+  side_sq = huge_scale * huge_scale;  //squared lenght of a side
+  diag_sq = 2.0 * side_sq;  //squared lenght of a diagonal
+  igl::squared_edge_lengths(V_huge,F,L_sq);
+  ASSERT_EQ(F.rows(), L_sq.rows());
+  ASSERT_EQ(3, L_sq.cols());
+  for(int f = 0;f<L_sq.rows();f++)
+  {
+    //All edge_lengths_squared must be exactly side_sq or diag_sq
+    for(int e = 0;e<3;e++)
+        if (L_sq(f,e) > 1.1*side_sq)
+           ASSERT_EQ(diag_sq, L_sq(f,e));
+        else
+           ASSERT_EQ(side_sq, L_sq(f,e));
+    //All sides sum exactly side_sq + side_sq + diag_sq
+    ASSERT_EQ(L_sq.row(f).sum(), side_sq + side_sq + diag_sq);
+  }
+ 
+  //Check the equilateral triangles
+  igl::squared_edge_lengths(V_huge,F_tet,L_sq);
+  ASSERT_EQ(F_tet.rows(), L_sq.rows());
+  ASSERT_EQ(3, L_sq.cols());
+  //All edges measure sqrt(2)
+  for(int f = 0;f<L_sq.rows();f++)
+  {
+    //All edge_lengths_squared must be exactly "diag_sq"
+    for(int e = 0;e<3;e++)
+       ASSERT_EQ(diag_sq, L_sq(f,e));
+  }
+
+  //Scale the cube to have tiny sides
+  side_sq = tiny_scale * tiny_scale;  //squared lenght of a side
+  diag_sq = 2.0 * side_sq;  //squared lenght of a diagonal
+  igl::squared_edge_lengths(V_tiny,F,L_sq);
+  ASSERT_EQ(F.rows(), L_sq.rows());
+  ASSERT_EQ(3, L_sq.cols());
+  for(int f = 0;f<L_sq.rows();f++)
+  {
+    //All edge_lengths_squared must be exactly side_sq or diag_sq
+    for(int e = 0;e<3;e++)
+        if (L_sq(f,e) > 1.1*side_sq)
+           ASSERT_EQ(diag_sq, L_sq(f,e));
+        else
+           ASSERT_EQ(side_sq, L_sq(f,e));
+    //All sides sum exactly side_sq + side_sq + diag_sq
+    ASSERT_EQ(L_sq.row(f).sum(), side_sq + side_sq + diag_sq);
+  }
+
+  //Check the regular tetrahedron
+  igl::squared_edge_lengths(V_tiny,F_tet,L_sq);
+  ASSERT_EQ(F_tet.rows(), L_sq.rows());
+  ASSERT_EQ(3, L_sq.cols());
+  //All edges measure sqrt(2)
+  for(int f = 0;f<L_sq.rows();f++)
+  {
+    //All edge_lengths_squared must be exactly "diag_sq"
+    for(int e = 0;e<3;e++)
+       ASSERT_EQ(diag_sq, L_sq(f,e));
+  }
+
+}

tests/include/igl/tet_tet_adjacency.cpp

@@ -0,0 +1,41 @@
+#include <test_common.h>
+#include <igl/tet_tet_adjacency.h>
+#include <igl/readMESH.h>
+#include <iostream>
+
+
+class tet_tet_adjacency : public ::testing::TestWithParam<std::string> {};
+
+TEST_P(tet_tet_adjacency, dot)
+{
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F, T, TT,TTi;
+  // Load example mesh: GetParam() will be name of mesh file
+  igl::readMESH(test_common::data_path(GetParam()), V, T, F);
+  igl::tet_tet_adjacency(T, TT, TTi);
+  ASSERT_EQ(T.rows(), TT.rows());
+  ASSERT_EQ(T.rows(), TTi.rows());
+  ASSERT_EQ(T.cols(),TT.cols());
+  ASSERT_EQ(T.cols(),TTi.cols());
+  for(int t = 0;t<T.rows();t++)
+  {
+    for(int c = 0; c<4 ;c++)
+    {
+      if(TT(t, c) >= 0)
+      {
+        ASSERT_LT(TT(t, c), T.rows());
+        ASSERT_GE(TTi(t, c), 0);
+        ASSERT_LT(TTi(t, c), 4);
+        ASSERT_EQ(TT(TT(t, c), TTi(t,c)) , t);
+      }
+    }
+  }
+}
+
+INSTANTIATE_TEST_CASE_P
+(
+  tet_meshes,
+  tet_tet_adjacency,
+  ::testing::ValuesIn(test_common::tet_meshes()),
+  test_common::string_test_name
+);

tests/include/igl/triangle_triangle_adjacency.cpp

@@ -0,0 +1,44 @@
+
+#include <test_common.h>
+#include <igl/triangle_triangle_adjacency.h>
+#include <Eigen/Geometry>
+
+class triangle_triangle_adjacency : public ::testing::TestWithParam<std::string> {};
+
+TEST_P(triangle_triangle_adjacency, dot)
+{
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F,TT,TTi;
+  // Load example mesh: GetParam() will be name of mesh file
+  test_common::load_mesh(GetParam(), V, F);
+  igl::triangle_triangle_adjacency(F,TT,TTi);
+  ASSERT_EQ(F.rows(),TT.rows());
+  ASSERT_EQ(F.rows(),TTi.rows());
+  ASSERT_EQ(F.cols(),TT.cols());
+  ASSERT_EQ(F.cols(),TTi.cols());
+  for(int f = 0;f<F.rows();f++)
+  {
+    for(int c = 0;c<3;c++)
+    {
+      if(TT(f,c) >= 0)
+      {
+        ASSERT_LT(TT(f,c) , F.rows());
+        ASSERT_GE(TTi(f,c) , 0);
+        ASSERT_LT(TTi(f,c) , 3);
+        ASSERT_EQ( TT(TT(f,c),TTi(f,c)) , f);
+      }
+    }
+  }
+  // ASSERT_EQ(a,b);
+  // ASSERT_TRUE(a==b);
+  // ASSERT_NEAR(a,b,1e-15)
+  // ASSERT_LT(a,1e-12);
+}
+
+INSTANTIATE_TEST_CASE_P
+(
+  manifold_meshes,
+  triangle_triangle_adjacency,
+  ::testing::ValuesIn(test_common::manifold_meshes()),
+  test_common::string_test_name
+);

tests/include/igl/unique.cpp

@@ -0,0 +1,84 @@
+#include <test_common.h>
+#include <igl/unique_rows.h>
+#include <igl/matrix_to_list.h>
+
+TEST(unique,matrix)
+{
+  Eigen::VectorXi A(12);
+  A = (Eigen::VectorXd::Random(A.size(),1).array().abs()*9).cast<int>();
+  Eigen::VectorXi C,IA,IC;
+  igl::unique_rows(A,C,IA,IC);
+  std::vector<bool> inA(A.maxCoeff()+1,false);
+  for(int i = 0;i<A.size();i++)
+  {
+    inA[A(i)] = true;
+    ASSERT_EQ(A(i),C(IC(i)));
+  }
+  std::vector<bool> inC(inA.size(),false);
+  // Expect a column vector
+  ASSERT_EQ(1,C.cols());
+  for(int i = 0;i<C.size();i++)
+  {
+    // Should be the first time finding this
+    ASSERT_FALSE(inC[C(i)]);
+    // Mark as found
+    inC[C(i)] = true;
+    // Should be something also found in A
+    ASSERT_TRUE(inA[C(i)]);
+    ASSERT_EQ(C(i),A(IA(i)));
+  }
+  for(int i = 0;i<inC.size();i++)
+  {
+    ASSERT_EQ(inC[i],inA[i]);
+  }
+}
+
+TEST(unique_rows,matrix)
+{
+  Eigen::MatrixXi A(50,4);
+  A = (Eigen::MatrixXi::Random(A.rows(),A.cols()).array().abs()*9).cast<int>();
+  Eigen::MatrixXi C;
+  Eigen::VectorXi IA,IC;
+  igl::unique_rows(A,C,IA,IC);
+  ASSERT_EQ(A.cols(),C.cols());
+  ASSERT_EQ(A.rows(),IC.size());
+  ASSERT_EQ(C.rows(),IA.size());
+  std::map<std::vector<int>,bool> inA;
+  for(int i = 0;i<A.rows();i++)
+  {
+    Eigen::RowVectorXi Ai = A.row(i);
+    std::vector<int> vAi;
+    igl::matrix_to_list(Ai,vAi);
+    inA[vAi] = true;
+    for(int j = 0;j<A.cols();j++)
+    {
+      ASSERT_EQ(A(i,j),C(IC(i),j));
+    }
+  }
+  std::map<std::vector<int>,bool> inC;
+  for(int i = 0;i<C.rows();i++)
+  {
+    Eigen::RowVectorXi Ci = C.row(i);
+    std::vector<int> vCi;
+    igl::matrix_to_list(Ci,vCi);
+    // Should be the first time finding this
+    ASSERT_FALSE(inC[vCi]);
+    // Mark as found
+    inC[vCi] = true;
+    // Should be something also found in A
+    ASSERT_TRUE(inA[vCi]);
+    for(int j = 0;j<A.cols();j++)
+    {
+      ASSERT_EQ(C(i,j),A(IA(i),j));
+    }
+  }
+  ASSERT_EQ(inC.size(),inA.size());
+  for(const auto pair : inA)
+  {
+    ASSERT_EQ(inC[pair.first],inA[pair.first]);
+  }
+  for(const auto pair : inC)
+  {
+    ASSERT_EQ(inC[pair.first],inA[pair.first]);
+  }
+}

tests/include/igl/upsample.cpp

@@ -0,0 +1,57 @@
+
+#include <test_common.h>
+#include <igl/upsample.h>
+
+class upsample : public ::testing::TestWithParam<std::string> {};
+
+TEST(upsample, single_triangle)
+{
+  Eigen::MatrixXi NF_groundtruth(4,3);
+  NF_groundtruth << 0,3,5 ,1,4,3 ,3,4,5 ,4,2,5;
+  Eigen::MatrixXd NV_groundtruth(6,2);
+  NV_groundtruth <<0,0 ,1,0 ,0,1 ,0.5,0 ,0.5,0.5 ,0,0.5;
+  Eigen::MatrixXd S_groundtruth(6,3);
+  S_groundtruth<<1,0,0 ,0,1,0 ,0,0,1 ,0.5,0.5,0 ,0,0.5,0.5 ,0.5,0,0.5;
+
+  Eigen::MatrixXi F(1,3);
+  F<<0,1,2;
+  Eigen::MatrixXd V(3,2);
+  V<<0,0,1,0,0,1;
+  Eigen::MatrixXi NF;
+  Eigen::MatrixXd NV;
+  Eigen::SparseMatrix<double> S;
+  igl::upsample(V.rows(),F,S,NF);
+  test_common::assert_eq(NF_groundtruth,NF);
+  test_common::assert_eq(S_groundtruth,Eigen::MatrixXd(S));
+  igl::upsample(V,F,NV,NF);
+  test_common::assert_eq(NF_groundtruth,NF);
+  test_common::assert_eq(NV_groundtruth,NV);
+}
+
+TEST_P(upsample, V_comes_first_F_ordering)
+{
+  Eigen::MatrixXd V,NV;
+  Eigen::MatrixXi F,NF;
+  // Load example mesh: GetParam() will be name of mesh file
+  test_common::load_mesh(GetParam(), V, F);
+  igl::upsample(V,F,NV,NF);
+  ASSERT_GE(NV.rows(),V.rows());
+  ASSERT_EQ(NF.rows(),4*F.rows());
+  // V should be first part of V
+  test_common::assert_eq(V,NV.topLeftCorner(V.rows(),V.cols()));
+  // Expect a particular order 
+  for(int f = 0;f<F.rows();f++)
+  {
+    ASSERT_EQ( F(f,0), NF((f*4)+0,0) );
+    ASSERT_EQ( F(f,1), NF((f*4)+1,0) );
+    ASSERT_EQ( F(f,2), NF((f*4)+3,1) );
+  }
+}
+
+INSTANTIATE_TEST_CASE_P
+(
+  manifold_meshes,
+  upsample,
+  ::testing::ValuesIn(test_common::manifold_meshes()),
+  test_common::string_test_name
+);

tests/test_common.h

@@ -0,0 +1,177 @@
+#pragma once
+
+
+#include <igl/read_triangle_mesh.h>
+#include <igl/find.h>
+#include <igl/readDMAT.h>
+
+#include <Eigen/Core>
+#include <gtest/gtest.h>
+
+#include <cctype>
+#include <string>
+#include <functional>
+#include <algorithm>
+#include <tuple>
+
+namespace test_common 
+{
+  // Input:
+  //   s  arbitrary string
+  // Returns s with all non-alphanumeric characters replaced with underscores '_'
+  inline std::string safe_test_name(std::string s)
+  {
+    std::for_each(s.begin(),s.end(),[](char &c){if(!std::isalnum(c)) c='_';});
+    return s;
+  };
+  inline std::string string_test_name(const ::testing::TestParamInfo<std::string>& info)
+  {
+    return test_common::safe_test_name(info.param);
+  };
+  inline std::vector<std::string> closed_genus_0_meshes()
+  {
+    return 
+    {
+      "cube.obj",
+      "decimated-knight.obj",
+      "boolean_minus_test_cube.obj",
+      "boolean_minus_test_green.obj",
+    };
+  };
+  inline std::vector<std::string> closed_manifold_meshes()
+  {
+    std::vector<std::string> meshes = closed_genus_0_meshes();
+    meshes.insert(meshes.end(),
+    {
+      "TinyTorus.obj",
+    });
+    return meshes;
+  };
+  inline std::vector<std::string> manifold_meshes()
+  {
+    std::vector<std::string> meshes = closed_manifold_meshes();
+    meshes.insert(meshes.end(),
+    {
+      "bunny.off",
+      "elephant.off",
+      "hemisphere.obj",
+    });
+    return meshes;
+  };
+  inline std::vector<std::string> tet_meshes()
+  {
+    return 
+    {
+      "decimated-knight.mesh"
+    };
+  };
+  inline std::vector<std::string> all_meshes()
+  {
+    std::vector<std::string> meshes = manifold_meshes();
+    meshes.insert(meshes.end(),
+    {
+      "truck.obj",
+    });
+    return meshes;
+  };
+  inline std::string data_path(std::string s)
+  {
+    return std::string(LIBIGL_DATA_DIR) + "/" + s;
+  };
+
+  // TODO: this seems like a pointless indirection. Should just find and
+  // replace test_common::load_mesh(X,...) with
+  // igl::read_triangle_mesh(test_common::data_path(X),...)
+  template<typename DerivedV, typename DerivedF>
+  void load_mesh(
+    const std::string& filename, 
+    Eigen::PlainObjectBase<DerivedV>& V,
+    Eigen::PlainObjectBase<DerivedF>& F)
+  {
+    igl::read_triangle_mesh(data_path(filename), V, F);
+  }
+
+  // TODO: this seems like a pointless indirection. Should just find and
+  // replace test_common::load_matrix(X,...) with
+  // igl::readDMAT(test_common::data_path(X),...)
+  template<typename Derived>
+  void load_matrix(
+    const std::string& filename,
+    Eigen::PlainObjectBase<Derived>& M) 
+  {
+    igl::readDMAT(data_path(filename), M);
+  }
+  template <typename DerivedA, typename DerivedB>
+  void assert_eq(
+    const Eigen::MatrixBase<DerivedA> & A,
+    const Eigen::MatrixBase<DerivedB> & B)
+  {
+    // Sizes should match
+    ASSERT_EQ(A.rows(),B.rows());
+    ASSERT_EQ(A.cols(),B.cols());
+    for(int i = 0;i<A.rows();i++)
+    {
+      for(int j = 0;j<A.cols();j++)
+      {
+        // Create an ijv tuple to trick GoogleTest into printing (i,j) so we
+        // know where the disagreement is.
+        std::tuple<int,int,typename DerivedA::Scalar> Aijv {i,j,A(i,j)};
+        std::tuple<int,int,typename DerivedB::Scalar> Bijv {i,j,B(i,j)};
+        ASSERT_EQ(Aijv,Bijv);
+      }
+    }
+  }
+  template <typename DerivedA, typename DerivedB>
+  void assert_eq(
+    const Eigen::SparseMatrix<DerivedA> & A,
+    const Eigen::SparseMatrix<DerivedB> & B)
+  {
+    // Sizes should match
+    ASSERT_EQ(A.rows(),B.rows());
+    ASSERT_EQ(A.cols(),B.cols());
+    Eigen::Matrix<long int,Eigen::Dynamic, 1> AI,AJ;
+    Eigen::Matrix<long int,Eigen::Dynamic, 1> BI,BJ;
+    Eigen::Matrix<DerivedA,Eigen::Dynamic, 1> AV;
+    Eigen::Matrix<DerivedB,Eigen::Dynamic, 1> BV;
+    // Assumes A and B are in same Major Ordering
+    igl::find(A,AI,AJ,AV);
+    igl::find(B,BI,BJ,BV);
+    // This doesn't generalized to assert_near nicely, and it makes it hard to
+    // tell which entries are different:
+    assert_eq(AI,BI);
+    assert_eq(AJ,BJ);
+    assert_eq(AV,BV);
+  }
+
+  template <typename DerivedA, typename DerivedB, typename EpsType>
+  void assert_near(
+    const Eigen::MatrixBase<DerivedA> & A,
+    const Eigen::MatrixBase<DerivedB> & B,
+    const EpsType & eps)
+  {
+    // Sizes should match
+    ASSERT_EQ(A.rows(),B.rows());
+    ASSERT_EQ(A.cols(),B.cols());
+    for(int i = 0;i<A.rows();i++)
+    {
+      for(int j = 0;j<A.cols();j++)
+      {
+        // Create an ijv tuple to trick GoogleTest into printing (i,j) so we
+        // know where the disagreement is.
+        //
+        // Equivalent to ASSERT_NEAR(Aijv,Bijv)
+        {
+          std::tuple<int,int,typename DerivedA::Scalar> Aijv {i,j,A(i,j)};
+          std::tuple<int,int,typename DerivedB::Scalar> Bijv {i,j,B(i,j)+eps};
+          ASSERT_LT(Aijv,Bijv);
+        }
+        {
+          std::tuple<int,int,typename DerivedA::Scalar> Aijv {i,j,A(i,j)+eps};
+          std::tuple<int,int,typename DerivedB::Scalar> Bijv {i,j,B(i,j)};
+          ASSERT_GT(Aijv,Bijv);
+        }
+      }
+    }
+  }
+
+}