diff --git a/cmake/igl/igl_add_tutorial.cmake b/cmake/igl/igl_add_tutorial.cmake
index e81513326..9077f4bb1 100644
--- a/cmake/igl/igl_add_tutorial.cmake
+++ b/cmake/igl/igl_add_tutorial.cmake
@@ -7,7 +7,9 @@ function(igl_add_tutorial name)
     endforeach()
 
     message(STATUS "Creating libigl tutorial: ${name}")
-    add_executable(${name} ${CMAKE_CURRENT_SOURCE_DIR}/${name}/main.cpp)
+    # get all cpp files in ${CMAKE_CURRENT_SOURCE_DIR}/${name}/
+    file(GLOB SRCFILES ${CMAKE_CURRENT_SOURCE_DIR}/${name}/*.cpp)
+    add_executable(${name} ${SRCFILES})
     target_link_libraries(${name} PRIVATE
         igl::core
         igl::tutorial_data
diff --git a/tutorial/406_FastAutomaticSkinningTransformations/main.cpp b/tutorial/406_FastAutomaticSkinningTransformations/main.cpp
index ad034248b..124912fa5 100644
--- a/tutorial/406_FastAutomaticSkinningTransformations/main.cpp
+++ b/tutorial/406_FastAutomaticSkinningTransformations/main.cpp
@@ -1,16 +1,6 @@
-#include <igl/colon.h>
-#include <igl/directed_edge_orientations.h>
-#include <igl/directed_edge_parents.h>
-#include <igl/forward_kinematics.h>
-#include <igl/PI.h>
-#include <igl/partition.h>
-#include <igl/max.h>
-#include <igl/lbs_matrix.h>
+#include "precomputation.h"
+
 #include <igl/slice.h>
-#include <igl/deform_skeleton.h>
-#include <igl/dqs.h>
-#include <igl/lbs_matrix.h>
-#include <igl/columnize.h>
 #include <igl/readDMAT.h>
 #include <igl/readOBJ.h>
 #include <igl/arap.h>
@@ -40,7 +30,6 @@ double bbd = 1.0;
 bool resolve = true;
 igl::ARAPData arap_data,arap_grouped_data;
 igl::ArapDOFData<Eigen::MatrixXd,double> arap_dof_data;
-Eigen::SparseMatrix<double> Aeq;
 
 enum ModeType
 {
@@ -149,57 +138,8 @@ int main(int argc, char *argv[])
   U=V;
   MatrixXd W;
   igl::readDMAT(TUTORIAL_SHARED_PATH "/armadillo-weights.dmat",W);
-  igl::lbs_matrix_column(V,W,M);
-
-  // Cluster according to weights
-  VectorXi G;
-  {
-    VectorXi S;
-    VectorXd D;
-    igl::partition(W,50,G,S,D);
-  }
-
-  // vertices corresponding to handles (those with maximum weight)
-  {
-    VectorXd maxW;
-    igl::max(W,1,maxW,b);
-  }
-
-  // Precomputation for FAST
-  cout<<"Initializing Fast Automatic Skinning Transformations..."<<endl;
-  // number of weights
-  const int m = W.cols();
-  Aeq.resize(m*3,m*3*(3+1));
-  vector<Triplet<double> > ijv;
-  for(int i = 0;i<m;i++)
-  {
-    RowVector4d homo;
-    homo << V.row(b(i)),1.;
-    for(int d = 0;d<3;d++)
-    {
-      for(int c = 0;c<(3+1);c++)
-      {
-        ijv.push_back(Triplet<double>(3*i + d,i + c*m*3 + d*m, homo(c)));
-      }
-    }
-  }
-  Aeq.setFromTriplets(ijv.begin(),ijv.end());
-  igl::arap_dof_precomputation(V,F,M,G,arap_dof_data);
-  igl::arap_dof_recomputation(VectorXi(),Aeq,arap_dof_data);
-  // Initialize
-  MatrixXd Istack = MatrixXd::Identity(3,3+1).replicate(1,m);
-  igl::columnize(Istack,m,2,L);
-
-  // Precomputation for ARAP
-  cout<<"Initializing ARAP..."<<endl;
-  arap_data.max_iter = 1;
-  igl::arap_precomputation(V,F,V.cols(),b,arap_data);
-  // Grouped arap
-  cout<<"Initializing ARAP with grouped edge-sets..."<<endl;
-  arap_grouped_data.max_iter = 2;
-  arap_grouped_data.G = G;
-  igl::arap_precomputation(V,F,V.cols(),b,arap_grouped_data);
 
+  precomputation(V,F,W,M,b,L,arap_data,arap_grouped_data,arap_dof_data);
 
   // bounding box diagonal
   bbd = (V.colwise().maxCoeff()- V.colwise().minCoeff()).norm();
diff --git a/tutorial/406_FastAutomaticSkinningTransformations/precomputation.cpp b/tutorial/406_FastAutomaticSkinningTransformations/precomputation.cpp
new file mode 100644
index 000000000..724b58380
--- /dev/null
+++ b/tutorial/406_FastAutomaticSkinningTransformations/precomputation.cpp
@@ -0,0 +1,76 @@
+#include "precomputation.h"
+
+#include <igl/arap.h>
+#include <igl/arap_dof.h>
+#include <igl/lbs_matrix.h>
+#include <igl/columnize.h>
+#include <igl/partition.h>
+#include <igl/max.h>
+
+
+void precomputation(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  const Eigen::MatrixXd & W,
+  Eigen::MatrixXd & M,
+  Eigen::VectorXi & b,
+  Eigen::MatrixXd & L,
+  igl::ARAPData & arap_data,
+  igl::ARAPData & arap_grouped_data,
+  igl::ArapDOFData<Eigen::MatrixXd,double> & arap_dof_data)
+{
+  using namespace Eigen;
+  using namespace std;
+  igl::lbs_matrix_column(V,W,M);
+
+  // Cluster according to weights
+  VectorXi G;
+  {
+    VectorXi S;
+    VectorXd D;
+    igl::partition(W,50,G,S,D);
+  }
+
+  // vertices corresponding to handles (those with maximum weight)
+  {
+    VectorXd maxW;
+    igl::max(W,1,maxW,b);
+  }
+
+  // Precomputation for FAST
+  cout<<"Initializing Fast Automatic Skinning Transformations..."<<endl;
+  // number of weights
+  const int m = W.cols();
+  Eigen::SparseMatrix<double> Aeq;
+  Aeq.resize(m*3,m*3*(3+1));
+  vector<Triplet<double> > ijv;
+  for(int i = 0;i<m;i++)
+  {
+    RowVector4d homo;
+    homo << V.row(b(i)),1.;
+    for(int d = 0;d<3;d++)
+    {
+      for(int c = 0;c<(3+1);c++)
+      {
+        ijv.push_back(Triplet<double>(3*i + d,i + c*m*3 + d*m, homo(c)));
+      }
+    }
+  }
+  Aeq.setFromTriplets(ijv.begin(),ijv.end());
+  igl::arap_dof_precomputation(V,F,M,G,arap_dof_data);
+  igl::arap_dof_recomputation(VectorXi(),Aeq,arap_dof_data);
+  // Initialize
+  MatrixXd Istack = MatrixXd::Identity(3,3+1).replicate(1,m);
+  igl::columnize(Istack,m,2,L);
+
+  // Precomputation for ARAP
+  cout<<"Initializing ARAP..."<<endl;
+  arap_data.max_iter = 1;
+  igl::arap_precomputation(V,F,V.cols(),b,arap_data);
+  // Grouped arap
+  cout<<"Initializing ARAP with grouped edge-sets..."<<endl;
+  arap_grouped_data.max_iter = 2;
+  arap_grouped_data.G = G;
+  igl::arap_precomputation(V,F,V.cols(),b,arap_grouped_data);
+
+}
diff --git a/tutorial/406_FastAutomaticSkinningTransformations/precomputation.h b/tutorial/406_FastAutomaticSkinningTransformations/precomputation.h
new file mode 100644
index 000000000..a31d2c41d
--- /dev/null
+++ b/tutorial/406_FastAutomaticSkinningTransformations/precomputation.h
@@ -0,0 +1,21 @@
+
+#include <Eigen/Core>
+
+namespace igl
+{
+  class ARAPData;
+  template <typename LbsMatrixType, typename SSCALAR>
+  // Seriously, windows? Why are struct and class different?
+  struct ArapDOFData;
+}
+
+void precomputation(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  const Eigen::MatrixXd & W,
+  Eigen::MatrixXd & M,
+  Eigen::VectorXi & b,
+  Eigen::MatrixXd & L,
+  igl::ARAPData & arap_data,
+  igl::ARAPData & arap_grouped_data,
+  igl::ArapDOFData<Eigen::MatrixXd,double> & arap_dof_data);
diff --git a/tutorial/610_CSGTree/get_mesh.cpp b/tutorial/610_CSGTree/get_mesh.cpp
new file mode 100644
index 000000000..dae9f4bae
--- /dev/null
+++ b/tutorial/610_CSGTree/get_mesh.cpp
@@ -0,0 +1,89 @@
+#include "get_mesh.h"
+#include <igl/copyleft/cgal/CSGTree.h>
+
+void get_mesh(
+    const Eigen::MatrixXd &VA,
+    const Eigen::MatrixXi &FA,
+    const Eigen::MatrixXd &VB,
+    const Eigen::MatrixXi &FB,
+    const Eigen::MatrixXd &VC,
+    const Eigen::MatrixXi &FC,
+    const Eigen::MatrixXd &VD,
+    const Eigen::MatrixXi &FD,
+    const Eigen::MatrixXd &VE,
+    const Eigen::MatrixXi &FE,
+    const int view_id,
+    Eigen::MatrixXd &V,
+    Eigen::MatrixXi &F,
+    Eigen::VectorXd &I)
+{
+  using namespace Eigen;
+  const auto & set_mesh = [&V,&F,&I](const Eigen::MatrixXd &V_, const Eigen::MatrixXi &F_, const int i)
+  {
+    V = V_;
+    F = F_;
+    I = VectorXd::Constant(F.rows(),1,i);
+  };
+  switch(view_id)
+  {
+    case 0:
+      set_mesh(VA,FA,5);
+      break;
+    case 1:
+      set_mesh(VB,FB,4);
+      break;
+    case 2:
+      set_mesh(VC,FC,3);
+      break;
+    case 3:
+      set_mesh(VD,FD,2);
+      break;
+    case 4:
+      set_mesh(VE,FE,1);
+      break;
+    default:
+    {
+      igl::copyleft::cgal::CSGTree M;
+      VectorXi J;
+      switch(view_id)
+      {
+        case 5:
+          // Compute result of (A ∩ B)
+          M = {{VA,FA},{VB,FB},"i"};
+          J = M.J().array()+0;
+          break;
+        case 6:
+          // Compute result of (C ∪ D)
+          M = {{VC,FC},{VD,FD},"u"};
+          J = M.J().array()+FA.rows()+FB.rows();
+          break;
+        case 7:
+          // Compute result of (C ∪ D) ∪ E
+          M = {{{VC,FC},{VD,FD},"u"},{VE,FE},"u"};
+          J = M.J().array()+FA.rows()+FB.rows();
+          break;
+        case 8:
+          // Compute result of (A ∩ B) \ ((C ∪ D) ∪ E)
+          M = {{{VA,FA},{VB,FB},"i"},{{{VC,FC},{VD,FD},"u"},{VE,FE},"u"},"m"};
+          J = M.J().array()+0;
+          break;
+        default:
+          assert(false && "unknown view id");
+      }
+      V = M.cast_V<MatrixXd>();
+      F = M.F();
+      I.resize(M.F().rows(),1);
+      // Compute colors based on original facets
+      for(int f = 0;f<M.F().rows();f++)
+      {
+        const int j = J(f);
+        I(f) = 
+          (int)(j<FA.rows())+
+          (int)(j<FA.rows()+FB.rows())+
+          (int)(j<FA.rows()+FB.rows()+FC.rows())+
+          (int)(j<FA.rows()+FB.rows()+FC.rows()+FD.rows())+
+          (int)(j<FA.rows()+FB.rows()+FC.rows()+FD.rows()+FE.rows());
+      }
+    }
+  }
+}
diff --git a/tutorial/610_CSGTree/get_mesh.h b/tutorial/610_CSGTree/get_mesh.h
new file mode 100644
index 000000000..934d9d67d
--- /dev/null
+++ b/tutorial/610_CSGTree/get_mesh.h
@@ -0,0 +1,18 @@
+#include <Eigen/Core>
+
+void get_mesh(
+    const Eigen::MatrixXd &VA,
+    const Eigen::MatrixXi &FA,
+    const Eigen::MatrixXd &VB,
+    const Eigen::MatrixXi &FB,
+    const Eigen::MatrixXd &VC,
+    const Eigen::MatrixXi &FC,
+    const Eigen::MatrixXd &VD,
+    const Eigen::MatrixXi &FD,
+    const Eigen::MatrixXd &VE,
+    const Eigen::MatrixXi &FE,
+    const int view_id,
+    Eigen::MatrixXd &V,
+    Eigen::MatrixXi &F,
+    Eigen::VectorXd &I);
+
diff --git a/tutorial/610_CSGTree/main.cpp b/tutorial/610_CSGTree/main.cpp
index a2e8cbdc6..e78b336ed 100644
--- a/tutorial/610_CSGTree/main.cpp
+++ b/tutorial/610_CSGTree/main.cpp
@@ -1,5 +1,5 @@
+#include "get_mesh.h"
 #include <igl/read_triangle_mesh.h>
-#include <igl/copyleft/cgal/CSGTree.h>
 #include <igl/opengl/glfw/Viewer.h>
 #include <igl/jet.h>
 #include <Eigen/Core>
@@ -8,7 +8,6 @@
 int main(int argc, char * argv[])
 {
   using namespace Eigen;
-  using namespace igl::copyleft::cgal;
   using namespace std;
   using namespace igl;
   cout<<R"(
@@ -30,76 +29,11 @@ int main(int argc, char * argv[])
   const auto & update = [&]()
   {
     viewer.data().clear();
-    // CSGTree templated on type of F
-    VectorXd I;
-    const auto & set_mesh = 
-      [&](const MatrixXd & V, const MatrixXi & F, const int i)
-    {
-      viewer.data().set_mesh(V,F);
-      I = VectorXd::Constant(F.rows(),1,i);
-    };
-    switch(view_id)
-    {
-      case 0:
-        set_mesh(VA,FA,5);
-        break;
-      case 1:
-        set_mesh(VB,FB,4);
-        break;
-      case 2:
-        set_mesh(VC,FC,3);
-        break;
-      case 3:
-        set_mesh(VD,FD,2);
-        break;
-      case 4:
-        set_mesh(VE,FE,1);
-        break;
-      default:
-      {
-        CSGTree M;
-        VectorXi J;
-        switch(view_id)
-        {
-          case 5:
-            // Compute result of (A ∩ B)
-            M = {{VA,FA},{VB,FB},"i"};
-            J = M.J().array()+0;
-            break;
-          case 6:
-            // Compute result of (C ∪ D)
-            M = {{VC,FC},{VD,FD},"u"};
-            J = M.J().array()+FA.rows()+FB.rows();
-            break;
-          case 7:
-            // Compute result of (C ∪ D) ∪ E
-            M = {{{VC,FC},{VD,FD},"u"},{VE,FE},"u"};
-            J = M.J().array()+FA.rows()+FB.rows();
-            break;
-          case 8:
-            // Compute result of (A ∩ B) \ ((C ∪ D) ∪ E)
-            M = {{{VA,FA},{VB,FB},"i"},{{{VC,FC},{VD,FD},"u"},{VE,FE},"u"},"m"};
-            J = M.J().array()+0;
-            break;
-          default:
-            assert(false && "unknown view id");
-        }
-        viewer.data().set_mesh(M.cast_V<MatrixXd>(),M.F());
-        I.resize(M.F().rows(),1);
-        // Compute colors based on original facets
-        for(int f = 0;f<M.F().rows();f++)
-        {
-          const int j = J(f);
-          I(f) = 
-            (int)(j<FA.rows())+
-            (int)(j<FA.rows()+FB.rows())+
-            (int)(j<FA.rows()+FB.rows()+FC.rows())+
-            (int)(j<FA.rows()+FB.rows()+FC.rows()+FD.rows())+
-            (int)(j<FA.rows()+FB.rows()+FC.rows()+FD.rows()+FE.rows());
-        }
-      }
-    }
-
+    Eigen::MatrixXd V;
+    Eigen::MatrixXi F;
+    Eigen::VectorXd I;
+    get_mesh(VA,FA,VB,FB,VC,FC,VD,FD,VE,FE,view_id,V,F,I);
+    viewer.data().set_mesh(V,F);
     MatrixXd C;
     jet(I,1,5,C);
     viewer.data().set_colors(C);
diff --git a/tutorial/709_SLIM/check_mesh_for_issues.cpp b/tutorial/709_SLIM/check_mesh_for_issues.cpp
new file mode 100644
index 000000000..9116b0e21
--- /dev/null
+++ b/tutorial/709_SLIM/check_mesh_for_issues.cpp
@@ -0,0 +1,44 @@
+#include "check_mesh_for_issues.h"
+#include <iostream>
+
+#include <igl/doublearea.h>
+#include <igl/vertex_components.h>
+#include <igl/euler_characteristic.h>
+#include <igl/is_edge_manifold.h>
+#include <igl/adjacency_matrix.h>
+#include <Eigen/Sparse>
+
+void check_mesh_for_issues(Eigen::MatrixXd& V, Eigen::MatrixXi& F) {
+
+  using namespace std;
+  Eigen::SparseMatrix<double> A;
+  igl::adjacency_matrix(F,A);
+
+  Eigen::MatrixXi C, Ci;
+  igl::vertex_components(A, C, Ci);
+
+  int connected_components = Ci.rows();
+  if (connected_components!=1) {
+    cout << "Error! Input has multiple connected components" << endl; exit(1);
+  }
+  int euler_char = igl::euler_characteristic(F);
+  if (euler_char!=1) 
+  {
+    cout << 
+      "Error! Input does not have a disk topology, it's euler char is " << 
+      euler_char << endl; 
+    exit(1);
+  }
+  bool is_edge_manifold = igl::is_edge_manifold(F);
+  if (!is_edge_manifold) {
+    cout << "Error! Input is not an edge manifold" << endl; exit(1);
+  }
+
+  Eigen::VectorXd areas; igl::doublearea(V,F,areas);
+  const double eps = 1e-14;
+  for (int i = 0; i < areas.rows(); i++) {
+    if (areas(i) < eps) {
+      cout << "Error! Input has zero area faces" << endl; exit(1);
+    }
+  }
+}
diff --git a/tutorial/709_SLIM/check_mesh_for_issues.h b/tutorial/709_SLIM/check_mesh_for_issues.h
new file mode 100644
index 000000000..116b184fb
--- /dev/null
+++ b/tutorial/709_SLIM/check_mesh_for_issues.h
@@ -0,0 +1,3 @@
+#include <Eigen/Core>
+
+void check_mesh_for_issues(Eigen::MatrixXd& V, Eigen::MatrixXi& F);
diff --git a/tutorial/709_SLIM/get_cube_corner_constraints.cpp b/tutorial/709_SLIM/get_cube_corner_constraints.cpp
new file mode 100644
index 000000000..71f29ca6a
--- /dev/null
+++ b/tutorial/709_SLIM/get_cube_corner_constraints.cpp
@@ -0,0 +1,90 @@
+
+#include "get_cube_corner_constraints.h"
+#include <igl/PI.h>
+#include <igl/cat.h>
+#include <set>
+
+void get_cube_corner_constraints(Eigen::MatrixXd& V, Eigen::MatrixXi& F, Eigen::VectorXi& b, Eigen::MatrixXd& bc) {
+  using namespace std;
+  double min_x,max_x,min_y,max_y,min_z,max_z;
+  min_x = V.col(0).minCoeff(); max_x = V.col(0).maxCoeff();
+  min_y = V.col(1).minCoeff(); max_y = V.col(1).maxCoeff();
+  min_z = V.col(2).minCoeff(); max_z = V.col(2).maxCoeff();
+
+
+  // get all cube corners
+  b.resize(8,1); bc.resize(8,3);
+  int x;
+  for (int i = 0; i < V.rows(); i++) {
+    if (V.row(i) == Eigen::RowVector3d(min_x,min_y,min_z)) b(0) = i;
+    if (V.row(i) == Eigen::RowVector3d(min_x,min_y,max_z)) b(1) = i;
+    if (V.row(i) == Eigen::RowVector3d(min_x,max_y,min_z)) b(2) = i;
+    if (V.row(i) == Eigen::RowVector3d(min_x,max_y,max_z)) b(3) = i;
+    if (V.row(i) == Eigen::RowVector3d(max_x,min_y,min_z)) b(4) = i;
+    if (V.row(i) == Eigen::RowVector3d(max_x,max_y,min_z)) b(5) = i;
+    if (V.row(i) == Eigen::RowVector3d(max_x,min_y,max_z)) b(6) = i;
+    if (V.row(i) == Eigen::RowVector3d(max_x,max_y,max_z)) b(7) = i;
+  }
+
+  // get all cube edges
+  std::set<int> cube_edge1; Eigen::VectorXi cube_edge1_vec;
+  for (int i = 0; i < V.rows(); i++) {
+    if ((V(i,0) == min_x && V(i,1) == min_y)) {
+      cube_edge1.insert(i);
+    }
+  }
+  Eigen::VectorXi edge1;
+  int_set_to_eigen_vector(cube_edge1, edge1);
+
+  std::set<int> cube_edge2; Eigen::VectorXi edge2;
+  for (int i = 0; i < V.rows(); i++) {
+    if ((V(i,0) == max_x && V(i,1) == max_y)) {
+      cube_edge2.insert(i);
+    }
+  }
+  int_set_to_eigen_vector(cube_edge2, edge2);
+  b = igl::cat(1,edge1,edge2);
+
+  std::set<int> cube_edge3; Eigen::VectorXi edge3;
+  for (int i = 0; i < V.rows(); i++) {
+    if ((V(i,0) == max_x && V(i,1) == min_y)) {
+      cube_edge3.insert(i);
+    }
+  }
+  int_set_to_eigen_vector(cube_edge3, edge3);
+  b = igl::cat(1,b,edge3);
+
+  std::set<int> cube_edge4; Eigen::VectorXi edge4;
+  for (int i = 0; i < V.rows(); i++) {
+    if ((V(i,0) == min_x && V(i,1) == max_y)) {
+      cube_edge4.insert(i);
+    }
+  }
+  int_set_to_eigen_vector(cube_edge4, edge4);
+  b = igl::cat(1,b,edge4);
+
+  bc.resize(b.rows(),3);
+  Eigen::Matrix3d m; m = Eigen::AngleAxisd(0.3 * igl::PI, Eigen::Vector3d(1./sqrt(2.),1./sqrt(2.),0.)/*Eigen::Vector3d::UnitX()*/);
+  int i = 0;
+  for (; i < cube_edge1.size(); i++) {
+    Eigen::RowVector3d edge_rot_center(min_x,min_y,(min_z+max_z)/2.);
+    bc.row(i) = (V.row(b(i)) - edge_rot_center) * m + edge_rot_center;
+  }
+  for (; i < cube_edge1.size() + cube_edge2.size(); i++) {
+    Eigen::RowVector3d edge_rot_center(max_x,max_y,(min_z+max_z)/2.);
+    bc.row(i) = (V.row(b(i)) - edge_rot_center) * m.transpose() + edge_rot_center;
+  }
+  for (; i < cube_edge1.size() + cube_edge2.size() + cube_edge3.size(); i++) {
+    bc.row(i) = 0.75*V.row(b(i));
+  }
+  for (; i < b.rows(); i++) {
+    bc.row(i) = 0.75*V.row(b(i));
+  }
+}
+
+void int_set_to_eigen_vector(const std::set<int>& int_set, Eigen::VectorXi& vec) {
+  vec.resize(int_set.size()); int idx = 0;
+  for(auto f : int_set) {
+      vec(idx) = f; idx++;
+    }
+}
diff --git a/tutorial/709_SLIM/get_cube_corner_constraints.h b/tutorial/709_SLIM/get_cube_corner_constraints.h
new file mode 100644
index 000000000..d55c2e564
--- /dev/null
+++ b/tutorial/709_SLIM/get_cube_corner_constraints.h
@@ -0,0 +1,5 @@
+
+#include <Eigen/Core>
+#include <set>
+void get_cube_corner_constraints(Eigen::MatrixXd& V_o, Eigen::MatrixXi& F, Eigen::VectorXi& b, Eigen::MatrixXd& bc);
+void int_set_to_eigen_vector(const std::set<int>& int_set, Eigen::VectorXi& vec);
diff --git a/tutorial/709_SLIM/get_soft_constraint_for_circle.cpp b/tutorial/709_SLIM/get_soft_constraint_for_circle.cpp
new file mode 100644
index 000000000..1be43b66e
--- /dev/null
+++ b/tutorial/709_SLIM/get_soft_constraint_for_circle.cpp
@@ -0,0 +1,32 @@
+#include "get_soft_constraint_for_circle.h"
+#include <igl/boundary_loop.h>
+
+void get_soft_constraint_for_circle(Eigen::MatrixXd& V_o, Eigen::MatrixXi& F, Eigen::VectorXi& b, Eigen::MatrixXd& bc) {
+
+    Eigen::VectorXi bnd;
+    igl::boundary_loop(F,bnd);
+    const int B_STEPS = 22; // constraint every B_STEPS vertices of the boundary
+
+    b.resize(bnd.rows()/B_STEPS);
+    bc.resize(b.rows(),2);
+
+    int c_idx = 0;
+    for (int i = B_STEPS; i < bnd.rows(); i+=B_STEPS) {
+        b(c_idx) = bnd(i);
+        c_idx++;
+    }
+
+    bc.row(0) = V_o.row(b(0)); // keep it there for now
+    bc.row(1) = V_o.row(b(2));
+    bc.row(2) = V_o.row(b(3));
+    bc.row(3) = V_o.row(b(4));
+    bc.row(4) = V_o.row(b(5));
+
+
+    bc.row(0) << V_o(b(0),0), 0;
+    bc.row(4) << V_o(b(4),0), 0;
+    bc.row(2) << V_o(b(2),0), 0.1;
+    bc.row(3) << V_o(b(3),0), 0.05;
+    bc.row(1) << V_o(b(1),0), -0.15;
+    bc.row(5) << V_o(b(5),0), +0.15;
+}
diff --git a/tutorial/709_SLIM/get_soft_constraint_for_circle.h b/tutorial/709_SLIM/get_soft_constraint_for_circle.h
new file mode 100644
index 000000000..9410057fb
--- /dev/null
+++ b/tutorial/709_SLIM/get_soft_constraint_for_circle.h
@@ -0,0 +1,3 @@
+#include <Eigen/Core>
+
+void get_soft_constraint_for_circle(Eigen::MatrixXd& V_o, Eigen::MatrixXi& F, Eigen::VectorXi& b, Eigen::MatrixXd& bc);
diff --git a/tutorial/709_SLIM/main.cpp b/tutorial/709_SLIM/main.cpp
index 0e43627d5..15e338d24 100644
--- a/tutorial/709_SLIM/main.cpp
+++ b/tutorial/709_SLIM/main.cpp
@@ -1,44 +1,29 @@
-#include <iostream>
+#include "check_mesh_for_issues.h"
+#include "get_soft_constraint_for_circle.h"
+#include "get_cube_corner_constraints.h"
+#include "param_2d_demo_iter.h"
+
 
 #include <igl/slim.h>
 
-#include <igl/vertex_components.h>
 #include <igl/readOBJ.h>
 #include <igl/writeOBJ.h>
 #include <igl/Timer.h>
 
-#include <igl/boundary_loop.h>
-#include <igl/map_vertices_to_circle.h>
-#include <igl/harmonic.h>
 #include <igl/MappingEnergyType.h>
 #include <igl/serialize.h>
 #include <igl/read_triangle_mesh.h>
 #include <igl/opengl/glfw/Viewer.h>
-#include <igl/flipped_triangles.h>
-#include <igl/euler_characteristic.h>
 #include <igl/barycenter.h>
-#include <igl/adjacency_matrix.h>
-#include <igl/is_edge_manifold.h>
-#include <igl/doublearea.h>
-#include <igl/cat.h>
-#include <igl/PI.h>
-
-#include <stdlib.h>
 
+#include <iostream>
+#include <cstdlib>
 #include <string>
 #include <vector>
 
-using namespace std;
-using namespace Eigen;
-
-void check_mesh_for_issues(Eigen::MatrixXd& V, Eigen::MatrixXi& F);
-void param_2d_demo_iter(igl::opengl::glfw::Viewer& viewer);
-void get_soft_constraint_for_circle(Eigen::MatrixXd& V_o, Eigen::MatrixXi& F, Eigen::VectorXi& b, Eigen::MatrixXd& bc);
 void soft_const_demo_iter(igl::opengl::glfw::Viewer& viewer);
 void deform_3d_demo_iter(igl::opengl::glfw::Viewer& viewer);
-void get_cube_corner_constraints(Eigen::MatrixXd& V_o, Eigen::MatrixXi& F, Eigen::VectorXi& b, Eigen::MatrixXd& bc);
 void display_3d_mesh(igl::opengl::glfw::Viewer& viewer);
-void int_set_to_eigen_vector(const std::set<int>& int_set, Eigen::VectorXi& vec);
 
 Eigen::MatrixXd V;
 Eigen::MatrixXi F;
@@ -59,7 +44,7 @@ bool key_down(igl::opengl::glfw::Viewer& viewer, unsigned char key, int modifier
   if (key == ' ') {
     switch (demo_type) {
       case PARAM_2D: {
-        param_2d_demo_iter(viewer);
+        param_2d_demo_iter(V,F,uv_scale_param,first_iter,sData,timer,viewer);
         break;
       }
       case SOFT_CONST: {
@@ -78,47 +63,8 @@ bool key_down(igl::opengl::glfw::Viewer& viewer, unsigned char key, int modifier
   return false;
 }
 
-void param_2d_demo_iter(igl::opengl::glfw::Viewer& viewer) {
-  if (first_iter) {
-    timer.start();
-    igl::read_triangle_mesh(TUTORIAL_SHARED_PATH "/face.obj", V, F);
-    check_mesh_for_issues(V,F);
-    cout << "\tMesh is valid!" << endl;
-
-    Eigen::MatrixXd uv_init;
-    Eigen::VectorXi bnd; Eigen::MatrixXd bnd_uv;
-    igl::boundary_loop(F,bnd);
-    igl::map_vertices_to_circle(V,bnd,bnd_uv);
-
-    igl::harmonic(V,F,bnd,bnd_uv,1,uv_init);
-    if (igl::flipped_triangles(uv_init,F).size() != 0) {
-      igl::harmonic(F,bnd,bnd_uv,1,uv_init); // use uniform laplacian
-    }
-
-    cout << "initialized parametrization" << endl;
-
-    sData.slim_energy = igl::MappingEnergyType::SYMMETRIC_DIRICHLET;
-    Eigen::VectorXi b; Eigen::MatrixXd bc;
-    slim_precompute(V,F,uv_init,sData, igl::MappingEnergyType::SYMMETRIC_DIRICHLET, b,bc,0);
-
-    uv_scale_param = 15 * (1./sqrt(sData.mesh_area));
-    viewer.data().set_mesh(V, F);
-    viewer.core().align_camera_center(V,F);
-    viewer.data().set_uv(sData.V_o*uv_scale_param);
-    viewer.data().compute_normals();
-    viewer.data().show_texture = true;
-
-    first_iter = false;
-  } else {
-    timer.start();
-    slim_solve(sData,1); // 1 iter
-    viewer.data().set_uv(sData.V_o*uv_scale_param);
-  }
-  cout << "time = " << timer.getElapsedTime() << endl;
-  cout << "energy = " << sData.energy << endl;
-}
-
 void soft_const_demo_iter(igl::opengl::glfw::Viewer& viewer) {
+  using namespace std;
   if (first_iter) {
 
     igl::read_triangle_mesh(TUTORIAL_SHARED_PATH "/circle.obj", V, F);
@@ -146,6 +92,7 @@ void soft_const_demo_iter(igl::opengl::glfw::Viewer& viewer) {
 }
 
 void deform_3d_demo_iter(igl::opengl::glfw::Viewer& viewer) {
+  using namespace std;
   if (first_iter) {
     timer.start();
     igl::readOBJ(TUTORIAL_SHARED_PATH "/cube_40k.obj", V, F);
@@ -172,6 +119,8 @@ void deform_3d_demo_iter(igl::opengl::glfw::Viewer& viewer) {
 }
 
 void display_3d_mesh(igl::opengl::glfw::Viewer& viewer) {
+  using namespace std;
+  using namespace Eigen;
   MatrixXd V_temp; MatrixXi F_temp;
   Eigen::MatrixXd Barycenters;
 
@@ -210,6 +159,8 @@ void display_3d_mesh(igl::opengl::glfw::Viewer& viewer) {
 }
 
 int main(int argc, char *argv[]) {
+  using namespace std;
+  using namespace Eigen;
 
   cerr << "Press space for running an iteration." << std::endl;
   cerr << "Syntax: " << argv[0] << " demo_number (1 to 3)" << std::endl;
@@ -257,150 +208,5 @@ int main(int argc, char *argv[]) {
   return 0;
 }
 
-void check_mesh_for_issues(Eigen::MatrixXd& V, Eigen::MatrixXi& F) {
-
-  Eigen::SparseMatrix<double> A;
-  igl::adjacency_matrix(F,A);
-
-  Eigen::MatrixXi C, Ci;
-  igl::vertex_components(A, C, Ci);
-
-  int connected_components = Ci.rows();
-  if (connected_components!=1) {
-    cout << "Error! Input has multiple connected components" << endl; exit(1);
-  }
-  int euler_char = igl::euler_characteristic(F);
-  if (euler_char!=1) 
-  {
-    cout << 
-      "Error! Input does not have a disk topology, it's euler char is " << 
-      euler_char << endl; 
-    exit(1);
-  }
-  bool is_edge_manifold = igl::is_edge_manifold(F);
-  if (!is_edge_manifold) {
-    cout << "Error! Input is not an edge manifold" << endl; exit(1);
-  }
-
-  Eigen::VectorXd areas; igl::doublearea(V,F,areas);
-  const double eps = 1e-14;
-  for (int i = 0; i < areas.rows(); i++) {
-    if (areas(i) < eps) {
-      cout << "Error! Input has zero area faces" << endl; exit(1);
-    }
-  }
-}
-
-void get_soft_constraint_for_circle(Eigen::MatrixXd& V_o, Eigen::MatrixXi& F, Eigen::VectorXi& b, Eigen::MatrixXd& bc) {
-
-    Eigen::VectorXi bnd;
-    igl::boundary_loop(F,bnd);
-    const int B_STEPS = 22; // constraint every B_STEPS vertices of the boundary
-
-    b.resize(bnd.rows()/B_STEPS);
-    bc.resize(b.rows(),2);
-
-    int c_idx = 0;
-    for (int i = B_STEPS; i < bnd.rows(); i+=B_STEPS) {
-        b(c_idx) = bnd(i);
-        c_idx++;
-    }
-
-    bc.row(0) = V_o.row(b(0)); // keep it there for now
-    bc.row(1) = V_o.row(b(2));
-    bc.row(2) = V_o.row(b(3));
-    bc.row(3) = V_o.row(b(4));
-    bc.row(4) = V_o.row(b(5));
 
 
-    bc.row(0) << V_o(b(0),0), 0;
-    bc.row(4) << V_o(b(4),0), 0;
-    bc.row(2) << V_o(b(2),0), 0.1;
-    bc.row(3) << V_o(b(3),0), 0.05;
-    bc.row(1) << V_o(b(1),0), -0.15;
-    bc.row(5) << V_o(b(5),0), +0.15;
-}
-
-void get_cube_corner_constraints(Eigen::MatrixXd& V, Eigen::MatrixXi& F, Eigen::VectorXi& b, Eigen::MatrixXd& bc) {
-  double min_x,max_x,min_y,max_y,min_z,max_z;
-  min_x = V.col(0).minCoeff(); max_x = V.col(0).maxCoeff();
-  min_y = V.col(1).minCoeff(); max_y = V.col(1).maxCoeff();
-  min_z = V.col(2).minCoeff(); max_z = V.col(2).maxCoeff();
-
-
-  // get all cube corners
-  b.resize(8,1); bc.resize(8,3);
-  int x;
-  for (int i = 0; i < V.rows(); i++) {
-    if (V.row(i) == Eigen::RowVector3d(min_x,min_y,min_z)) b(0) = i;
-    if (V.row(i) == Eigen::RowVector3d(min_x,min_y,max_z)) b(1) = i;
-    if (V.row(i) == Eigen::RowVector3d(min_x,max_y,min_z)) b(2) = i;
-    if (V.row(i) == Eigen::RowVector3d(min_x,max_y,max_z)) b(3) = i;
-    if (V.row(i) == Eigen::RowVector3d(max_x,min_y,min_z)) b(4) = i;
-    if (V.row(i) == Eigen::RowVector3d(max_x,max_y,min_z)) b(5) = i;
-    if (V.row(i) == Eigen::RowVector3d(max_x,min_y,max_z)) b(6) = i;
-    if (V.row(i) == Eigen::RowVector3d(max_x,max_y,max_z)) b(7) = i;
-  }
-
-  // get all cube edges
-  std::set<int> cube_edge1; Eigen::VectorXi cube_edge1_vec;
-  for (int i = 0; i < V.rows(); i++) {
-    if ((V(i,0) == min_x && V(i,1) == min_y)) {
-      cube_edge1.insert(i);
-    }
-  }
-  Eigen::VectorXi edge1;
-  int_set_to_eigen_vector(cube_edge1, edge1);
-
-  std::set<int> cube_edge2; Eigen::VectorXi edge2;
-  for (int i = 0; i < V.rows(); i++) {
-    if ((V(i,0) == max_x && V(i,1) == max_y)) {
-      cube_edge2.insert(i);
-    }
-  }
-  int_set_to_eigen_vector(cube_edge2, edge2);
-  b = igl::cat(1,edge1,edge2);
-
-  std::set<int> cube_edge3; Eigen::VectorXi edge3;
-  for (int i = 0; i < V.rows(); i++) {
-    if ((V(i,0) == max_x && V(i,1) == min_y)) {
-      cube_edge3.insert(i);
-    }
-  }
-  int_set_to_eigen_vector(cube_edge3, edge3);
-  b = igl::cat(1,b,edge3);
-
-  std::set<int> cube_edge4; Eigen::VectorXi edge4;
-  for (int i = 0; i < V.rows(); i++) {
-    if ((V(i,0) == min_x && V(i,1) == max_y)) {
-      cube_edge4.insert(i);
-    }
-  }
-  int_set_to_eigen_vector(cube_edge4, edge4);
-  b = igl::cat(1,b,edge4);
-
-  bc.resize(b.rows(),3);
-  Eigen::Matrix3d m; m = Eigen::AngleAxisd(0.3 * igl::PI, Eigen::Vector3d(1./sqrt(2.),1./sqrt(2.),0.)/*Eigen::Vector3d::UnitX()*/);
-  int i = 0;
-  for (; i < cube_edge1.size(); i++) {
-    Eigen::RowVector3d edge_rot_center(min_x,min_y,(min_z+max_z)/2.);
-    bc.row(i) = (V.row(b(i)) - edge_rot_center) * m + edge_rot_center;
-  }
-  for (; i < cube_edge1.size() + cube_edge2.size(); i++) {
-    Eigen::RowVector3d edge_rot_center(max_x,max_y,(min_z+max_z)/2.);
-    bc.row(i) = (V.row(b(i)) - edge_rot_center) * m.transpose() + edge_rot_center;
-  }
-  for (; i < cube_edge1.size() + cube_edge2.size() + cube_edge3.size(); i++) {
-    bc.row(i) = 0.75*V.row(b(i));
-  }
-  for (; i < b.rows(); i++) {
-    bc.row(i) = 0.75*V.row(b(i));
-  }
-}
-
-void int_set_to_eigen_vector(const std::set<int>& int_set, Eigen::VectorXi& vec) {
-  vec.resize(int_set.size()); int idx = 0;
-  for(auto f : int_set) {
-      vec(idx) = f; idx++;
-    }
-}
diff --git a/tutorial/709_SLIM/param_2d_demo_iter.cpp b/tutorial/709_SLIM/param_2d_demo_iter.cpp
new file mode 100644
index 000000000..e8ed2429f
--- /dev/null
+++ b/tutorial/709_SLIM/param_2d_demo_iter.cpp
@@ -0,0 +1,61 @@
+#include "param_2d_demo_iter.h"
+#include "check_mesh_for_issues.h"
+#include <igl/read_triangle_mesh.h>
+#include <igl/opengl/glfw/Viewer.h>
+#include <igl/flipped_triangles.h>
+#include <igl/MappingEnergyType.h>
+#include <igl/harmonic.h>
+#include <igl/map_vertices_to_circle.h>
+#include <igl/boundary_loop.h>
+#include <igl/slim.h>
+#include <igl/Timer.h>
+
+void param_2d_demo_iter(
+    Eigen::MatrixXd & V,
+    Eigen::MatrixXi & F,
+    double & uv_scale_param,
+    bool & first_iter,
+    igl::SLIMData& sData,
+    igl::Timer & timer,
+    igl::opengl::glfw::Viewer& viewer)
+{
+  using namespace std;
+  if (first_iter) {
+    timer.start();
+    igl::read_triangle_mesh(TUTORIAL_SHARED_PATH "/face.obj", V, F);
+    check_mesh_for_issues(V,F);
+    cout << "\tMesh is valid!" << endl;
+
+    Eigen::MatrixXd uv_init;
+    Eigen::VectorXi bnd; Eigen::MatrixXd bnd_uv;
+    igl::boundary_loop(F,bnd);
+    igl::map_vertices_to_circle(V,bnd,bnd_uv);
+
+    igl::harmonic(V,F,bnd,bnd_uv,1,uv_init);
+    if (igl::flipped_triangles(uv_init,F).size() != 0) {
+      igl::harmonic(F,bnd,bnd_uv,1,uv_init); // use uniform laplacian
+    }
+
+    cout << "initialized parametrization" << endl;
+
+    sData.slim_energy = igl::MappingEnergyType::SYMMETRIC_DIRICHLET;
+    Eigen::VectorXi b; Eigen::MatrixXd bc;
+    slim_precompute(V,F,uv_init,sData, igl::MappingEnergyType::SYMMETRIC_DIRICHLET, b,bc,0);
+
+    uv_scale_param = 15 * (1./sqrt(sData.mesh_area));
+    viewer.data().set_mesh(V, F);
+    viewer.core().align_camera_center(V,F);
+    viewer.data().set_uv(sData.V_o*uv_scale_param);
+    viewer.data().compute_normals();
+    viewer.data().show_texture = true;
+
+    first_iter = false;
+  } else {
+    timer.start();
+    slim_solve(sData,1); // 1 iter
+    viewer.data().set_uv(sData.V_o*uv_scale_param);
+  }
+  cout << "time = " << timer.getElapsedTime() << endl;
+  cout << "energy = " << sData.energy << endl;
+}
+
diff --git a/tutorial/709_SLIM/param_2d_demo_iter.h b/tutorial/709_SLIM/param_2d_demo_iter.h
new file mode 100644
index 000000000..40915a505
--- /dev/null
+++ b/tutorial/709_SLIM/param_2d_demo_iter.h
@@ -0,0 +1,20 @@
+#include <Eigen/Core>
+namespace igl {
+  class SLIMData;
+  class Timer;
+  namespace opengl {
+    namespace glfw {
+      class Viewer;
+    }
+  }
+}
+
+void param_2d_demo_iter(
+    Eigen::MatrixXd & V,
+    Eigen::MatrixXi & F,
+    double & uv_scale_param,
+    bool & first_iter,
+    igl::SLIMData& sData,
+    igl::Timer & timer,
+    igl::opengl::glfw::Viewer& viewer);
+
diff --git a/tutorial/805_MeshImplicitFunction/contours.cpp b/tutorial/805_MeshImplicitFunction/contours.cpp
new file mode 100644
index 000000000..a0674723f
--- /dev/null
+++ b/tutorial/805_MeshImplicitFunction/contours.cpp
@@ -0,0 +1,202 @@
+#include "contours.h"
+#include <igl/unique_simplices.h>
+#include <igl/dual_contouring.h>
+#include <igl/get_seconds.h>
+#include <igl/grid.h>
+#include <igl/marching_cubes.h>
+#include <igl/per_corner_normals.h>
+#include <igl/parallel_for.h>
+#include <igl/per_corner_normals.h>
+#include <igl/per_face_normals.h>
+#include <igl/polygon_corners.h>
+#include <igl/slice.h>
+#include <igl/sparse_voxel_grid.h>
+
+void contours(
+    Eigen::MatrixXd & V,
+    Eigen::MatrixXi & E,
+    Eigen::MatrixXd & VV,
+    Eigen::MatrixXi & FF,
+    Eigen::MatrixXd & NN,
+    Eigen::MatrixXd & mcV,
+    Eigen::MatrixXi & mcF,
+    Eigen::MatrixXd & mcN)
+{
+  const auto & tictoc = []()
+  {
+    static double t_start = igl::get_seconds();
+    double diff = igl::get_seconds()-t_start;
+    t_start += diff;
+    return diff;
+  };
+
+  // Create an interesting shape with sharp features using SDF CSG with spheres.
+  const auto & sphere = [](
+      const Eigen::RowVector3d & c,
+      const double r,
+      const Eigen::RowVector3d & x)->double
+  {
+    return (x-c).norm() - r;
+  };
+  const std::function<double(const Eigen::RowVector3d & x)> f = 
+    [&](const Eigen::RowVector3d & x)->double
+  {
+    return 
+      std::min(
+      std::min(std::max(
+      sphere(Eigen::RowVector3d(-0.2,0,-0.2),0.5,x),
+      -sphere(Eigen::RowVector3d(+0.2,0,0.2),0.5,x)),
+       sphere(Eigen::RowVector3d(-0.15,0,-0.15),0.3,x)
+      ),
+      std::max(
+      std::max(
+        sphere(Eigen::RowVector3d(-0.2,-0.5,-0.2),0.6,x),x(1)+0.45),-0.6-x(1))
+      );
+  };
+  Eigen::RowVector3d p0(-0.2,0.5,-0.2);
+  assert(abs(f(p0)) < 1e-10 && "p0 should be on zero level-set");
+
+  // Simple finite difference gradients
+  const auto & fd = [](
+    const std::function<double(const Eigen::RowVector3d&)> & f,
+    const Eigen::RowVector3d & x)
+  {
+    const double eps = 1e-10;
+    Eigen::RowVector3d g;
+    for(int c = 0;c<3;c++)
+    {
+      const Eigen::RowVector3d xp = x+eps*Eigen::RowVector3d(c==0,c==1,c==2);
+      const double fp = f(xp);
+      const Eigen::RowVector3d xn = x-eps*Eigen::RowVector3d(c==0,c==1,c==2);
+      const double fn = f(xn);
+      g(c) = (fp-fn)/(2*eps);
+    }
+    return g;
+  };
+  const auto & f_grad = [&fd,&f](const Eigen::RowVector3d & x)
+  {
+    return fd(f,x).normalized();
+  };
+
+  // Quads
+  Eigen::MatrixXi Q;
+
+  // Grid parameters
+  const Eigen::RowVector3d min_corner(-2,-2,-2);
+  const Eigen::RowVector3d max_corner(+2,+2,+2);
+  const int s = 256;
+  int nx = s+1;
+  int ny = s+1;
+  int nz = s+1;
+  const Eigen::RowVector3d step =
+    (max_corner-min_corner).array()/(Eigen::RowVector3d(nx,ny,nz).array()-1);
+  // Sparse grid below assumes regular grid
+  assert((step(0) == step(1))&&(step(0) == step(2)));
+
+  // Dual contouring parameters
+  bool constrained = false;
+  bool triangles = false;
+  bool root_finding = true;
+  for(int pass = 0;pass<2;pass++)
+  {
+    const bool sparse = pass == 1;
+    printf("Using %s grid..\n",sparse?"sparse":"dense");
+    if(sparse)
+    {
+      // igl::sparse_voxel_grid assumes (0,0,0) lies on the grid. But dense igl::grid
+      // below won't necessarily do that depending on nx,ny,nz.
+      tictoc();
+      Eigen::MatrixXd GV;
+      Eigen::VectorXd Gf;
+      Eigen::Matrix<int,Eigen::Dynamic,8> GI;
+      igl::sparse_voxel_grid(p0,f,step(0),16.*pow(step(0),-2.),Gf,GV,GI);
+      const auto t_Gf = tictoc();
+      printf("  %5f secs to populate sparse grid of %ld cells\n",t_Gf+tictoc(),GI.rows());
+      // Dual contouring requires list of sparse edges (not cells)
+      // extract _all_ edges from sparse_voxel_grid (conservative)
+      Eigen::Matrix<int,Eigen::Dynamic,2> GI2;
+      {
+        Eigen::Matrix<int,Eigen::Dynamic,2> all_GI2(GI.rows()*12,2);
+        all_GI2 << 
+          // front
+          GI.col(0),GI.col(1),
+          GI.col(1),GI.col(2),
+          GI.col(2),GI.col(3),
+          GI.col(3),GI.col(0),
+          // back
+          GI.col(4+0),GI.col(4+1),
+          GI.col(4+1),GI.col(4+2),
+          GI.col(4+2),GI.col(4+3),
+          GI.col(4+3),GI.col(4+0),
+          // sides
+          GI.col(0),GI.col(4+0), 
+          GI.col(1),GI.col(4+1), 
+          GI.col(2),GI.col(4+2),
+          GI.col(3),GI.col(4+3);
+        Eigen::VectorXi _1,_2;
+        igl::unique_simplices(all_GI2,GI2,_1,_2);
+      }
+      tictoc();
+      Eigen::RowVector3d step =
+        (max_corner-min_corner).array()/(Eigen::RowVector3d(nx,ny,nz).array()-1);
+      igl::dual_contouring(
+        f,f_grad,step,Gf,GV,GI2,constrained,triangles,root_finding,V,Q);
+      printf("  %5f secs dual contouring\n",t_Gf+tictoc());
+      tictoc();
+      igl::marching_cubes(Gf,GV,GI,0.0,mcV,mcF);
+      printf("  %5f secs marching cubes\n",t_Gf+tictoc());
+    }else
+    {
+
+      tictoc();
+      igl::dual_contouring(
+        f,f_grad,min_corner,max_corner,nx,ny,nz,constrained,triangles,root_finding,V,Q);
+      printf("  %5f secs dual contouring\n",tictoc());
+      // build and sample grid
+      tictoc();
+      Eigen::MatrixXd GV;
+      igl::grid(Eigen::RowVector3i(nx,ny,nz),GV);
+      Eigen::VectorXd Gf(GV.rows());
+      igl::parallel_for(GV.rows(),[&](const int i)
+      {
+        GV.row(i).array() *= (max_corner-min_corner).array();
+        GV.row(i) += min_corner;
+        Gf(i) = f(GV.row(i));
+      },1000ul);
+      const auto t_grid = tictoc();
+      igl::marching_cubes(Gf,GV,nx,ny,nz,0,mcV,mcF);
+      const auto t_mc = tictoc();
+      printf("  %5f secs (%5f + %5f) marching cubes\n",t_grid+t_mc,t_grid,t_mc);
+    }
+  }
+
+  // Crisp (as possible) rendering of resulting MC triangle mesh
+  igl::per_corner_normals(mcV,mcF,20,mcN);
+  // Crisp rendering of resulting DC quad mesh with edges
+  Eigen::MatrixXd N;
+  Eigen::VectorXi J;
+  if(triangles)
+  {
+    VV = V;
+    FF = Q;
+    E.resize(Q.rows()*3,2);
+    E<<
+      Q.col(0), Q.col(1), 
+      Q.col(1), Q.col(2), 
+      Q.col(2), Q.col(0);
+  }else
+  {
+    Eigen::VectorXi I,C;
+    igl::polygon_corners(Q,I,C);
+    E.resize(Q.rows()*4,2);
+    E<<
+      Q.col(0), Q.col(1), 
+      Q.col(1), Q.col(2), 
+      Q.col(2), Q.col(3), 
+      Q.col(3), Q.col(0);
+    igl::per_face_normals(V,I,C,N,VV,FF,J);
+    igl::slice(N,J,1,NN);
+    igl::per_corner_normals(V,I,C,20,N,VV,FF,J,NN);
+  }
+
+}
diff --git a/tutorial/805_MeshImplicitFunction/contours.h b/tutorial/805_MeshImplicitFunction/contours.h
new file mode 100644
index 000000000..ab9d7d89b
--- /dev/null
+++ b/tutorial/805_MeshImplicitFunction/contours.h
@@ -0,0 +1,11 @@
+#include <Eigen/Core>
+
+void contours(
+  Eigen::MatrixXd & V,
+  Eigen::MatrixXi & E,
+  Eigen::MatrixXd & VV,
+  Eigen::MatrixXi & FF,
+  Eigen::MatrixXd & NN,
+  Eigen::MatrixXd & mcV,
+  Eigen::MatrixXi & mcF,
+  Eigen::MatrixXd & mcN);
diff --git a/tutorial/805_MeshImplicitFunction/main.cpp b/tutorial/805_MeshImplicitFunction/main.cpp
index 83f4b0667..e8a94ae5e 100644
--- a/tutorial/805_MeshImplicitFunction/main.cpp
+++ b/tutorial/805_MeshImplicitFunction/main.cpp
@@ -1,199 +1,17 @@
-#include <igl/unique_simplices.h>
-#include <igl/dual_contouring.h>
-#include <igl/get_seconds.h>
-#include <igl/grid.h>
-#include <igl/marching_cubes.h>
-#include <igl/per_corner_normals.h>
-#include <igl/parallel_for.h>
-#include <igl/per_corner_normals.h>
-#include <igl/per_face_normals.h>
-#include <igl/polygon_corners.h>
-#include <igl/slice.h>
-#include <igl/sparse_voxel_grid.h>
+#include "contours.h"
 #include <igl/opengl/glfw/Viewer.h>
 
 int main(int argc, char * argv[])
 {
-  const auto & tictoc = []()
-  {
-    static double t_start = igl::get_seconds();
-    double diff = igl::get_seconds()-t_start;
-    t_start += diff;
-    return diff;
-  };
-
-  // Create an interesting shape with sharp features using SDF CSG with spheres.
-  const auto & sphere = [](
-      const Eigen::RowVector3d & c,
-      const double r,
-      const Eigen::RowVector3d & x)->double
-  {
-    return (x-c).norm() - r;
-  };
-  const std::function<double(const Eigen::RowVector3d & x)> f = 
-    [&](const Eigen::RowVector3d & x)->double
-  {
-    return 
-      std::min(
-      std::min(std::max(
-      sphere(Eigen::RowVector3d(-0.2,0,-0.2),0.5,x),
-      -sphere(Eigen::RowVector3d(+0.2,0,0.2),0.5,x)),
-       sphere(Eigen::RowVector3d(-0.15,0,-0.15),0.3,x)
-      ),
-      std::max(
-      std::max(
-        sphere(Eigen::RowVector3d(-0.2,-0.5,-0.2),0.6,x),x(1)+0.45),-0.6-x(1))
-      );
-  };
-  Eigen::RowVector3d p0(-0.2,0.5,-0.2);
-  assert(abs(f(p0)) < 1e-10 && "p0 should be on zero level-set");
-
-  // Simple finite difference gradients
-  const auto & fd = [](
-    const std::function<double(const Eigen::RowVector3d&)> & f,
-    const Eigen::RowVector3d & x)
-  {
-    const double eps = 1e-10;
-    Eigen::RowVector3d g;
-    for(int c = 0;c<3;c++)
-    {
-      const Eigen::RowVector3d xp = x+eps*Eigen::RowVector3d(c==0,c==1,c==2);
-      const double fp = f(xp);
-      const Eigen::RowVector3d xn = x-eps*Eigen::RowVector3d(c==0,c==1,c==2);
-      const double fn = f(xn);
-      g(c) = (fp-fn)/(2*eps);
-    }
-    return g;
-  };
-  const auto & f_grad = [&fd,&f](const Eigen::RowVector3d & x)
-  {
-    return fd(f,x).normalized();
-  };
-
   Eigen::MatrixXd V;
-  Eigen::MatrixXi Q,F;
-  Eigen::MatrixXd mcV,mcN;
-  Eigen::MatrixXi mcF;
-
-  // Grid parameters
-  const Eigen::RowVector3d min_corner(-2,-2,-2);
-  const Eigen::RowVector3d max_corner(+2,+2,+2);
-  const int s = 256;
-  int nx = s+1;
-  int ny = s+1;
-  int nz = s+1;
-  const Eigen::RowVector3d step =
-    (max_corner-min_corner).array()/(Eigen::RowVector3d(nx,ny,nz).array()-1);
-  // Sparse grid below assumes regular grid
-  assert((step(0) == step(1))&&(step(0) == step(2)));
-
-  // Dual contouring parameters
-  bool constrained = false;
-  bool triangles = false;
-  bool root_finding = true;
-  for(int pass = 0;pass<2;pass++)
-  {
-    const bool sparse = pass == 1;
-    printf("Using %s grid..\n",sparse?"sparse":"dense");
-    if(sparse)
-    {
-      // igl::sparse_voxel_grid assumes (0,0,0) lies on the grid. But dense igl::grid
-      // below won't necessarily do that depending on nx,ny,nz.
-      tictoc();
-      Eigen::MatrixXd GV;
-      Eigen::VectorXd Gf;
-      Eigen::Matrix<int,Eigen::Dynamic,8> GI;
-      igl::sparse_voxel_grid(p0,f,step(0),16.*pow(step(0),-2.),Gf,GV,GI);
-      const auto t_Gf = tictoc();
-      printf("  %5f secs to populate sparse grid of %ld cells\n",t_Gf+tictoc(),GI.rows());
-      // Dual contouring requires list of sparse edges (not cells)
-      // extract _all_ edges from sparse_voxel_grid (conservative)
-      Eigen::Matrix<int,Eigen::Dynamic,2> GI2;
-      {
-        Eigen::Matrix<int,Eigen::Dynamic,2> all_GI2(GI.rows()*12,2);
-        all_GI2 << 
-          // front
-          GI.col(0),GI.col(1),
-          GI.col(1),GI.col(2),
-          GI.col(2),GI.col(3),
-          GI.col(3),GI.col(0),
-          // back
-          GI.col(4+0),GI.col(4+1),
-          GI.col(4+1),GI.col(4+2),
-          GI.col(4+2),GI.col(4+3),
-          GI.col(4+3),GI.col(4+0),
-          // sides
-          GI.col(0),GI.col(4+0), 
-          GI.col(1),GI.col(4+1), 
-          GI.col(2),GI.col(4+2),
-          GI.col(3),GI.col(4+3);
-        Eigen::VectorXi _1,_2;
-        igl::unique_simplices(all_GI2,GI2,_1,_2);
-      }
-      tictoc();
-      Eigen::RowVector3d step =
-        (max_corner-min_corner).array()/(Eigen::RowVector3d(nx,ny,nz).array()-1);
-      igl::dual_contouring(
-        f,f_grad,step,Gf,GV,GI2,constrained,triangles,root_finding,V,Q);
-      printf("  %5f secs dual contouring\n",t_Gf+tictoc());
-      tictoc();
-      igl::marching_cubes(Gf,GV,GI,0.0,mcV,mcF);
-      printf("  %5f secs marching cubes\n",t_Gf+tictoc());
-    }else
-    {
-
-      tictoc();
-      igl::dual_contouring(
-        f,f_grad,min_corner,max_corner,nx,ny,nz,constrained,triangles,root_finding,V,Q);
-      printf("  %5f secs dual contouring\n",tictoc());
-      // build and sample grid
-      tictoc();
-      Eigen::MatrixXd GV;
-      igl::grid(Eigen::RowVector3i(nx,ny,nz),GV);
-      Eigen::VectorXd Gf(GV.rows());
-      igl::parallel_for(GV.rows(),[&](const int i)
-      {
-        GV.row(i).array() *= (max_corner-min_corner).array();
-        GV.row(i) += min_corner;
-        Gf(i) = f(GV.row(i));
-      },1000ul);
-      const auto t_grid = tictoc();
-      igl::marching_cubes(Gf,GV,nx,ny,nz,0,mcV,mcF);
-      const auto t_mc = tictoc();
-      printf("  %5f secs (%5f + %5f) marching cubes\n",t_grid+t_mc,t_grid,t_mc);
-    }
-  }
-
-  // Crisp (as possible) rendering of resulting MC triangle mesh
-  igl::per_corner_normals(mcV,mcF,20,mcN);
-  // Crisp rendering of resulting DC quad mesh with edges
   Eigen::MatrixXi E;
-  Eigen::MatrixXd VV,N,NN;
-  Eigen::VectorXi J;
+  Eigen::MatrixXd VV;
   Eigen::MatrixXi FF;
-  if(triangles)
-  {
-    VV = V;
-    FF = Q;
-    E.resize(Q.rows()*3,2);
-    E<<
-      Q.col(0), Q.col(1), 
-      Q.col(1), Q.col(2), 
-      Q.col(2), Q.col(0);
-  }else
-  {
-    Eigen::VectorXi I,C;
-    igl::polygon_corners(Q,I,C);
-    E.resize(Q.rows()*4,2);
-    E<<
-      Q.col(0), Q.col(1), 
-      Q.col(1), Q.col(2), 
-      Q.col(2), Q.col(3), 
-      Q.col(3), Q.col(0);
-    igl::per_face_normals(V,I,C,N,VV,FF,J);
-    igl::slice(N,J,1,NN);
-    igl::per_corner_normals(V,I,C,20,N,VV,FF,J,NN);
-  }
+  Eigen::MatrixXd NN;
+  Eigen::MatrixXd mcV;
+  Eigen::MatrixXi mcF;
+  Eigen::MatrixXd mcN;
+  contours(V,E,VV,FF,NN,mcV,mcF,mcN);
 
   igl::opengl::glfw::Viewer vr;
   bool show_edges = true;
diff --git a/tutorial/806_HeatGeodesics/isolines_colormap.cpp b/tutorial/806_HeatGeodesics/isolines_colormap.cpp
new file mode 100644
index 000000000..dbeff64bf
--- /dev/null
+++ b/tutorial/806_HeatGeodesics/isolines_colormap.cpp
@@ -0,0 +1,19 @@
+#include "isolines_colormap.h"
+#include <igl/isolines_map.h>
+
+
+Eigen::MatrixXd isolines_colormap()
+{
+  const int num_intervals = 30;
+  Eigen::MatrixXd CM(num_intervals,3);
+  // Colormap texture
+  for(int i = 0;i<num_intervals;i++)
+  {
+    double t = double(num_intervals - i - 1)/double(num_intervals-1);
+    CM(i,0) = std::max(std::min(2.0*t-0.0,1.0),0.0);
+    CM(i,1) = std::max(std::min(2.0*t-1.0,1.0),0.0);
+    CM(i,2) = std::max(std::min(6.0*t-5.0,1.0),0.0);
+  }
+  igl::isolines_map(Eigen::MatrixXd(CM),CM);
+  return CM;
+}
diff --git a/tutorial/806_HeatGeodesics/isolines_colormap.h b/tutorial/806_HeatGeodesics/isolines_colormap.h
new file mode 100644
index 000000000..a55746150
--- /dev/null
+++ b/tutorial/806_HeatGeodesics/isolines_colormap.h
@@ -0,0 +1,3 @@
+#include <Eigen/Core>
+
+Eigen::MatrixXd isolines_colormap();
diff --git a/tutorial/806_HeatGeodesics/main.cpp b/tutorial/806_HeatGeodesics/main.cpp
index 554ef5a22..eb0a7c53e 100644
--- a/tutorial/806_HeatGeodesics/main.cpp
+++ b/tutorial/806_HeatGeodesics/main.cpp
@@ -1,40 +1,20 @@
+#include "isolines_colormap.h"
+#include "update.h"
 #include <igl/read_triangle_mesh.h>
-#include <igl/triangulated_grid.h>
 #include <igl/heat_geodesics.h>
-#include <igl/unproject_onto_mesh.h>
 #include <igl/avg_edge_length.h>
-#include <igl/isolines_map.h>
 #include <igl/opengl/glfw/Viewer.h>
-#include <igl/opengl/create_shader_program.h>
-#include <igl/opengl/destroy_shader_program.h>
 #include <iostream>
 
-void set_colormap(igl::opengl::glfw::Viewer & viewer)
-{
-  const int num_intervals = 30;
-  Eigen::MatrixXd CM(num_intervals,3);
-  // Colormap texture
-  for(int i = 0;i<num_intervals;i++)
-  {
-    double t = double(num_intervals - i - 1)/double(num_intervals-1);
-    CM(i,0) = std::max(std::min(2.0*t-0.0,1.0),0.0);
-    CM(i,1) = std::max(std::min(2.0*t-1.0,1.0),0.0);
-    CM(i,2) = std::max(std::min(6.0*t-5.0,1.0),0.0);
-  }
-  igl::isolines_map(Eigen::MatrixXd(CM),CM);
-  viewer.data().set_colormap(CM);
-}
-
 int main(int argc, char *argv[])
 {
-  // Create the peak height field
   Eigen::MatrixXi F;
   Eigen::MatrixXd V;
   igl::read_triangle_mesh( argc>1?argv[1]: TUTORIAL_SHARED_PATH "/beetle.off",V,F);
+  double t = std::pow(igl::avg_edge_length(V,F),2);
 
   // Precomputation
   igl::HeatGeodesicsData<double> data;
-  double t = std::pow(igl::avg_edge_length(V,F),2);
   const auto precompute = [&]()
   {
     if(!igl::heat_geodesics_precompute(V,F,t,data))
@@ -49,40 +29,15 @@ int main(int argc, char *argv[])
   bool down_on_mesh = false;
   const auto update = [&]()->bool
   {
-    int fid;
-    Eigen::Vector3f bc;
-    // Cast a ray in the view direction starting from the mouse position
-    double x = viewer.current_mouse_x;
-    double y = viewer.core().viewport(3) - viewer.current_mouse_y;
-    if(igl::unproject_onto_mesh(Eigen::Vector2f(x,y), viewer.core().view,
-      viewer.core().proj, viewer.core().viewport, V, F, fid, bc))
+    const double x = viewer.current_mouse_x;
+    const double y = viewer.core().viewport(3) - viewer.current_mouse_y;
+    Eigen::VectorXd D;
+    if(::update(
+      V,F,t,x,y,
+      viewer.core().view,viewer.core().proj,viewer.core().viewport,
+      data,
+      D))
     {
-      Eigen::VectorXd D;
-      // if big mesh, just use closest vertex. Otherwise, blend distances to
-      // vertices of face using barycentric coordinates.
-      if(F.rows()>100000)
-      {
-        // 3d position of hit
-        const Eigen::RowVector3d m3 =
-          V.row(F(fid,0))*bc(0) + V.row(F(fid,1))*bc(1) + V.row(F(fid,2))*bc(2);
-        int cid = 0;
-        Eigen::Vector3d(
-            (V.row(F(fid,0))-m3).squaredNorm(),
-            (V.row(F(fid,1))-m3).squaredNorm(),
-            (V.row(F(fid,2))-m3).squaredNorm()).minCoeff(&cid);
-        const int vid = F(fid,cid);
-        igl::heat_geodesics_solve(data,(Eigen::VectorXi(1,1)<<vid).finished(),D);
-      }else
-      {
-        D = Eigen::VectorXd::Zero(V.rows());
-        for(int cid = 0;cid<3;cid++)
-        {
-          const int vid = F(fid,cid);
-          Eigen::VectorXd Dc;
-          igl::heat_geodesics_solve(data,(Eigen::VectorXi(1,1)<<vid).finished(),Dc);
-          D += Dc*bc(cid);
-        }
-      }
       viewer.data().set_data(D);
       return true;
     }
@@ -150,7 +105,7 @@ int main(int argc, char *argv[])
   // Show mesh
   viewer.data().set_mesh(V, F);
   viewer.data().set_data(Eigen::VectorXd::Zero(V.rows()));
-  set_colormap(viewer);
+  viewer.data().set_colormap(isolines_colormap());
   viewer.data().show_lines = false;
   viewer.launch();
 
diff --git a/tutorial/806_HeatGeodesics/update.cpp b/tutorial/806_HeatGeodesics/update.cpp
new file mode 100644
index 000000000..ea1d1826a
--- /dev/null
+++ b/tutorial/806_HeatGeodesics/update.cpp
@@ -0,0 +1,51 @@
+#include "update.h"
+#include <igl/heat_geodesics.h>
+#include <igl/unproject_onto_mesh.h>
+
+bool update(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  const double t,
+  const double x,
+  const double y,
+  const Eigen::Matrix4f& view,
+  const Eigen::Matrix4f& proj,
+  const Eigen::Vector4f& viewport,
+  const igl::HeatGeodesicsData<double>& data,
+  Eigen::VectorXd& D)
+{
+  int fid;
+  Eigen::Vector3f bc;
+  // Cast a ray in the view direction starting from the mouse position
+  if(igl::unproject_onto_mesh(Eigen::Vector2f(x,y), view,
+    proj, viewport, V, F, fid, bc))
+  {
+    // if big mesh, just use closest vertex. Otherwise, blend distances to
+    // vertices of face using barycentric coordinates.
+    if(F.rows()>100000)
+    {
+      // 3d position of hit
+      const Eigen::RowVector3d m3 =
+        V.row(F(fid,0))*bc(0) + V.row(F(fid,1))*bc(1) + V.row(F(fid,2))*bc(2);
+      int cid = 0;
+      Eigen::Vector3d(
+          (V.row(F(fid,0))-m3).squaredNorm(),
+          (V.row(F(fid,1))-m3).squaredNorm(),
+          (V.row(F(fid,2))-m3).squaredNorm()).minCoeff(&cid);
+      const int vid = F(fid,cid);
+      igl::heat_geodesics_solve(data,(Eigen::VectorXi(1,1)<<vid).finished(),D);
+    }else
+    {
+      D = Eigen::VectorXd::Zero(V.rows());
+      for(int cid = 0;cid<3;cid++)
+      {
+        const int vid = F(fid,cid);
+        Eigen::VectorXd Dc;
+        igl::heat_geodesics_solve(data,(Eigen::VectorXi(1,1)<<vid).finished(),Dc);
+        D += Dc*bc(cid);
+      }
+    }
+    return true;
+  }
+  return false;
+}
diff --git a/tutorial/806_HeatGeodesics/update.h b/tutorial/806_HeatGeodesics/update.h
new file mode 100644
index 000000000..9af581c56
--- /dev/null
+++ b/tutorial/806_HeatGeodesics/update.h
@@ -0,0 +1,18 @@
+#include <Eigen/Core>
+namespace igl
+{
+  template <typename Scalar>
+  struct HeatGeodesicsData;
+}
+
+bool update(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  const double t,
+  const double x,
+  const double y,
+  const Eigen::Matrix4f& model,
+  const Eigen::Matrix4f& proj,
+  const Eigen::Vector4f& viewport,
+  const igl::HeatGeodesicsData<double>& data,
+  Eigen::VectorXd& D);