laplace equation example

Former-commit-id: dc1978ab82efea11a3f69457fa96525888843898

include/igl/frame_field_deformer.cpp

@@ -284,7 +284,7 @@ IGL_INLINE void Frame_field_deformer::compute_optimal_positions()
   for (int i=0;i<nfree;i++)
   {
     b.row(i) << 0.0, 0.0, 0.0;
-    for (int k=0;k<VT[i].size();k++)					// for all incident triangles
+    for (int k=0;k<(int)VT[i].size();k++)					// for all incident triangles
     {
       t = VT[i][k];												// incident tri
 			vi = (i==F(t,0))?0:(i==F(t,1))?1:(i==F(t,2))?2:3;	// index of i in t
@@ -345,7 +345,7 @@ IGL_INLINE void Frame_field_deformer::compute_optimal_positions()
   using namespace Eigen;
 
   WW.resize(F.rows());
-	for (int i=0;i<FF.size();i++)
+	for (int i=0;i<(int)FF.size();i++)
 	{
 		Vector3d v0,v1,v2;
 		v0 = FF[i].col(0);
@@ -360,7 +360,7 @@ IGL_INLINE void Frame_field_deformer::compute_optimal_positions()
 
 		// polar decomposition to discard rotational component (unnecessary but makes it easier)
 		Eigen::JacobiSVD<Matrix<double,3,3> > svd(AI, Eigen::ComputeFullU | Eigen::ComputeFullV );
-		Matrix<double,3,3>  au = svd.matrixU();
+		//Matrix<double,3,3>  au = svd.matrixU();
 		Matrix<double,3,3>  av = svd.matrixV();
 		DiagonalMatrix<double,3>	as(svd.singularValues());
 		WW[i] = av*as*av.transpose();

include/igl/intersect.cpp

@@ -10,6 +10,7 @@ template <class M>
 IGL_INLINE void igl::intersect(const M & A, const M & B, M & C)
 {
   // Stupid O(size(A) * size(B)) to do it
+  // Alec: This should be implemented by using unique and sort like `setdiff`
   M dyn_C(A.size() > B.size() ? A.size() : B.size(),1);
   // count of intersects
   int c = 0;

include/igl/is_border_vertex.h

@@ -20,7 +20,7 @@ namespace igl
   // Inputs:
   //   V  #V by dim list of vertex positions 
   //   F  #F by 3 list of triangle indices
-  // Returns vector of indices of vertices on open boundary of F
+  // Returns #V vector of bools revealing whether vertices are on boundary
   //
   // Known Bugs: does not depend on V
   // 

include/igl/list_to_matrix.cpp

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

include/igl/matrix_to_list.cpp

@@ -35,9 +35,12 @@ IGL_INLINE void igl::matrix_to_list(
   using namespace std;
   V.resize(M.size());
   // loop over rows
-  for(int i = 0;i<M.size();i++)
+  for(int j = 0;j<M.cols();j++)
   {
-    V[i] = M[i];
+    for(int i = 0;i<M.rows();i++)
+    {
+      V[i+j*M.rows()] = M(i,j);
+    }
   }
 }
 
@@ -61,4 +64,6 @@ template std::vector<Eigen::Matrix<int, -1, 1, 0, -1, 1>::Scalar, std::allocator
 template void igl::matrix_to_list<Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, std::vector<std::vector<Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, std::allocator<Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar> >, std::allocator<std::vector<Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, std::allocator<Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar> > > >&);
 template void igl::matrix_to_list<Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<std::vector<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Scalar, std::allocator<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Scalar> >, std::allocator<std::vector<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Scalar, std::allocator<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Scalar> > > >&);
 template void igl::matrix_to_list<Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, std::vector<Eigen::Matrix<double, -1, 1, 0, -1, 1>::Scalar, std::allocator<Eigen::Matrix<double, -1, 1, 0, -1, 1>::Scalar> >&);
+template void igl::matrix_to_list<Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Scalar, std::allocator<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Scalar> >&);
+template void igl::matrix_to_list<Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, std::vector<Eigen::Matrix<int, -1, 1, 0, -1, 1>::Scalar, std::allocator<Eigen::Matrix<int, -1, 1, 0, -1, 1>::Scalar> >&);
 #endif

include/igl/sort.cpp

@@ -232,4 +232,5 @@ template void igl::sort<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<in
 template void igl::sort<Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::sort<Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 2, 0, -1, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&);
 template void igl::sort<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::sort<Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/unique.cpp

@@ -11,6 +11,7 @@
 #include "SortableRow.h"
 #include "sortrows.h"
 #include "list_to_matrix.h"
+#include "matrix_to_list.h"
 
 #include <algorithm>
 #include <iostream>
@@ -61,6 +62,29 @@ IGL_INLINE void igl::unique(
 
 }
 
+template <
+  typename DerivedA,
+  typename DerivedC,
+  typename DerivedIA,
+  typename DerivedIC>
+IGL_INLINE void igl::unique(
+    const Eigen::PlainObjectBase<DerivedA> & A,
+    Eigen::PlainObjectBase<DerivedC> & C,
+    Eigen::PlainObjectBase<DerivedIA> & IA,
+    Eigen::PlainObjectBase<DerivedIC> & IC)
+{
+  using namespace std;
+  using namespace Eigen;
+  vector<typename DerivedA::Scalar > vA;
+  vector<typename DerivedC::Scalar > vC;
+  vector<size_t> vIA,vIC;
+  matrix_to_list(A,vA);
+  unique(vA,vC,vIA,vIC);
+  list_to_matrix(vC,C);
+  list_to_matrix(vIA,IA);
+  list_to_matrix(vIC,IC);
+}
+
 // Obsolete slow version converting to vectors
 // template <typename DerivedA, typename DerivedIA, typename DerivedIC>
 // IGL_INLINE void igl::unique_rows(
@@ -196,10 +220,13 @@ IGL_INLINE void igl::unique_rows(
 }
 
 #ifndef IGL_HEADER_ONLY
+// Explicit template specialization
 template void igl::unique<int>(std::vector<int, std::allocator<int> > const&, std::vector<int, std::allocator<int> >&, std::vector<size_t, std::allocator<size_t> >&, std::vector<size_t, std::allocator<size_t> >&);
 template void igl::unique_rows<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<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> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::unique_rows<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template void igl::unique_rows<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<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> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::unique_rows<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::unique_rows<Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template void igl::unique<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<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> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template void igl::unique<Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<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> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/unique.h

@@ -29,7 +29,16 @@ namespace igl
     std::vector<T> & C,
     std::vector<size_t> & IA,
     std::vector<size_t> & IC);
-
+  template <
+    typename DerivedA,
+    typename DerivedC,
+    typename DerivedIA,
+    typename DerivedIC>
+  IGL_INLINE void unique(
+      const Eigen::PlainObjectBase<DerivedA> & A,
+      Eigen::PlainObjectBase<DerivedC> & C,
+      Eigen::PlainObjectBase<DerivedIA> & IA,
+      Eigen::PlainObjectBase<DerivedIC> & IC);
   // Act like matlab's [C,IA,IC] = unique(X,'rows')
   //
   // Templates:

tutorial/303_LaplaceEquation/CMakeLists.txt

@@ -0,0 +1,11 @@
+cmake_minimum_required(VERSION 2.6)
+project(303_LaplaceEquation)
+
+include("../CMakeLists.shared")
+
+set(SOURCES
+${PROJECT_SOURCE_DIR}/main.cpp
+)
+
+add_executable(${CMAKE_PROJECT_NAME} ${SOURCES} ${SHARED_SOURCES})
+target_link_libraries(${CMAKE_PROJECT_NAME} ${SHARED_LIBRARIES})

tutorial/303_LaplaceEquation/main.cpp

@@ -0,0 +1,71 @@
+#include <igl/boundary_faces.h>
+#include <igl/colon.h>
+#include <igl/cotmatrix.h>
+#include <igl/jet.h>
+#include <igl/min_quad_with_fixed.h>
+#include <igl/readOFF.h>
+#include <igl/setdiff.h>
+#include <igl/slice.h>
+#include <igl/slice_into.h>
+#include <igl/unique.h>
+#include <igl/viewer/Viewer.h>
+#include <Eigen/Sparse>
+#include <iostream>
+int main(int argc, char *argv[])
+{
+  using namespace Eigen;
+  using namespace std;
+  MatrixXd V;
+  MatrixXi F;
+  igl::readOFF("../shared/camelhead.off",V,F);
+  // Find boundary edges
+  MatrixXi E;
+  igl::boundary_faces(F,E);
+  // Find boundary vertices
+  VectorXi b,IA,IC;
+  igl::unique(E,b,IA,IC);
+  // List of all vertex indices
+  VectorXi all,in;
+  igl::colon<int>(0,V.rows()-1,all);
+  // List of interior indices
+  igl::setdiff(all,b,in,IA);
+
+  // Construct and slice up Laplacian
+  SparseMatrix<double> L,L_in_in,L_in_b;
+  igl::cotmatrix(V,F,L);
+  igl::slice(L,in,in,L_in_in);
+  igl::slice(L,in,b,L_in_b);
+
+  // Dirichelet boundary conditions from z-coordinate
+  VectorXd bc;
+  VectorXd Z = V.col(2);
+  igl::slice(Z,b,bc);
+
+  // Solve PDE
+  SimplicialLLT<SparseMatrix<double > > solver(-L_in_in);
+  VectorXd Z_in = solver.solve(L_in_b*bc);
+  // slice into solution
+  igl::slice_into(Z_in,in,Z);
+
+  // Alternative, short hand
+  igl::min_quad_with_fixed_data<double> mqwf;
+  // Linear term is 0
+  VectorXd B = VectorXd::Zero(V.rows(),1);
+  // Empty constraints
+  VectorXd Beq;
+  SparseMatrix<double> Aeq;
+  // Our cotmatrix is _negative_ definite, so flip sign
+  igl::min_quad_with_fixed_precompute((-L).eval(),b,Aeq,true,mqwf);
+  igl::min_quad_with_fixed_solve(mqwf,B,bc,Beq,Z);
+
+  // Pseudo-color based on solution
+  MatrixXd C;
+  igl::jet(Z,true,C);
+
+  // Plot the mesh with pseudocolors
+  igl::Viewer viewer;
+  viewer.set_mesh(V, F);
+  viewer.options.show_lines = false;
+  viewer.set_colors(C);
+  viewer.launch();
+}

tutorial/cmake/FindLIBIGL.cmake

@@ -48,6 +48,13 @@ if(LIBIGL_USE_STATIC_LIBRARY)
     set(LIBIGL_FOUND FALSE)
   endif(NOT LIBIGLMATLAB_LIBRARY)
   set(LIBIGL_LIBRARIES ${LIBIGL_LIBRARIES}  ${LIBIGLMATLAB_LIBRARY})
+  find_package(Matlab REQUIRED)
+  if(MATLAB_FOUND)
+    set(LIBIGL_INCLUDE_DIR ${LIBIGL_INCLUDE_DIR}  ${MATLAB_INCLUDE_DIR})
+    set(LIBIGL_LIBRARIES ${LIBIGL_LIBRARIES}  ${MATLAB_LIBRARIES})
+  else(MATLAB_FOUND)
+    set(LIBIGL_FOUND FALSE)
+  endif(MATLAB_FOUND)
   FIND_LIBRARY( LIBIGLSVD3X3_LIBRARY NAMES iglsvd3x3 PATHS ${LIBIGL_LIB_DIRS})
   if(NOT LIBIGLSVD3X3_LIBRARY)
     set(LIBIGL_FOUND FALSE)

tutorial/cmake/FindMATLAB.cmake

@@ -1,54 +1,210 @@
+# - this module looks for Matlab
+# Defines:
+#  MATLAB_INCLUDE_DIR: include path for mex.h, engine.h
+#  MATLAB_LIBRARIES:   required libraries: libmex, etc
+#  MATLAB_MEX_LIBRARY: path to libmex.lib
+#  MATLAB_MX_LIBRARY:  path to libmx.lib
+#  MATLAB_MAT_LIBRARY:  path to libmat.lib # added
+#  MATLAB_ENG_LIBRARY: path to libeng.lib
+#  MATLAB_ROOT: path to Matlab's root directory
+
+# This file is part of Gerardus
+#
+# This is a derivative work of file FindMatlab.cmake released with
+# CMake v2.8, because the original seems to be a bit outdated and
+# doesn't work with my Windows XP and Visual Studio 10 install
 #
-# Try to find a MATLAB installation
-# Once done this will define
+# (Note that the original file does work for Ubuntu Natty)
 #
-# MATLAB_FOUND
-# MATLAB_INCLUDE_DIR
-# MATLAB_LIBRARIES
+# Author: Ramon Casero <rcasero@gmail.com>, Tom Doel
+# Version: 0.2.3
+# $Rev$
+# $Date$
 #
+# The original file was copied from an Ubuntu Linux install
+# /usr/share/cmake-2.8/Modules/FindMatlab.cmake
 
-FIND_PATH(MATLAB_INCLUDE_DIR engine.h
-      PATHS
-      /Applications/MATLAB_R2012b.app/extern/include
-      /Applications/MATLAB_R2013a.app/extern/include
-      /Applications/MATLAB_R2013b.app/extern/include
-      /Applications/MATLAB_R2014a.app/extern/include
-      /Applications/MATLAB_R2014b.app/extern/include
-      NO_DEFAULT_PATH)
-
-FIND_LIBRARY(MATLAB_LIBRARY1 eng
-  PATHS
-    /Applications/MATLAB_R2012b.app/bin/maci64
-    /Applications/MATLAB_R2013a.app/bin/maci64
-    /Applications/MATLAB_R2013b.app/bin/maci64
-    /Applications/MATLAB_R2014a.app/bin/maci64
-    /Applications/MATLAB_R2014b.app/bin/maci64
-  PATH_SUFFIXES`
-    a
-    lib64
-    lib
-    dylib
-    NO_DEFAULT_PATH
-)
+#=============================================================================
+# Copyright 2005-2009 Kitware, Inc.
+#
+# Distributed under the OSI-approved BSD License (the "License");
+# see accompanying file Copyright.txt for details.
+#
+# This software is distributed WITHOUT ANY WARRANTY; without even the
+# implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+# See the License for more information.
+#=============================================================================
+# (To distribute this file outside of CMake, substitute the full
+#  License text for the above reference.)
+
+SET(MATLAB_FOUND 0)
+IF(WIN32)
+  # Search for a version of Matlab available, starting from the most modern one to older versions
+  FOREACH(MATVER "7.14" "7.11" "7.10" "7.9" "7.8" "7.7" "7.6" "7.5" "7.4")
+    IF((NOT DEFINED MATLAB_ROOT) 
+        OR ("${MATLAB_ROOT}" STREQUAL "")
+        OR ("${MATLAB_ROOT}" STREQUAL "/registry"))
+      GET_FILENAME_COMPONENT(MATLAB_ROOT
+        "[HKEY_LOCAL_MACHINE\\SOFTWARE\\MathWorks\\MATLAB\\${MATVER};MATLABROOT]"
+        ABSOLUTE)
+      SET(MATLAB_VERSION ${MATVER})
+    ENDIF((NOT DEFINED MATLAB_ROOT) 
+      OR ("${MATLAB_ROOT}" STREQUAL "")
+      OR ("${MATLAB_ROOT}" STREQUAL "/registry"))
+  ENDFOREACH(MATVER)
+  
+  # Directory name depending on whether the Windows architecture is 32
+  # bit or 64 bit
+  set(CMAKE_SIZEOF_VOID_P 8) # Note: For some weird reason this variable is undefined in my system...
+  IF(CMAKE_SIZEOF_VOID_P MATCHES "4")
+    SET(WINDIR "win32")
+  ELSEIF(CMAKE_SIZEOF_VOID_P MATCHES "8")
+    SET(WINDIR "win64")
+  ELSE(CMAKE_SIZEOF_VOID_P MATCHES "4")
+    MESSAGE(FATAL_ERROR 
+      "CMAKE_SIZEOF_VOID_P (${CMAKE_SIZEOF_VOID_P}) doesn't indicate a valid platform")
+  ENDIF(CMAKE_SIZEOF_VOID_P MATCHES "4")
+
+  # Folder where the MEX libraries are, depending of the Windows compiler
+  IF(${CMAKE_GENERATOR} MATCHES "Visual Studio 6")
+    SET(MATLAB_LIBRARIES_DIR "${MATLAB_ROOT}/extern/lib/${WINDIR}/microsoft/msvc60")
+  ELSEIF(${CMAKE_GENERATOR} MATCHES "Visual Studio 7")
+    # Assume people are generally using Visual Studio 7.1,
+    # if using 7.0 need to link to: ../extern/lib/${WINDIR}/microsoft/msvc70
+    SET(MATLAB_LIBRARIES_DIR "${MATLAB_ROOT}/extern/lib/${WINDIR}/microsoft/msvc71")
+    # SET(MATLAB_LIBRARIES_DIR "${MATLAB_ROOT}/extern/lib/${WINDIR}/microsoft/msvc70")
+  ELSEIF(${CMAKE_GENERATOR} MATCHES "Borland")
+    # Assume people are generally using Borland 5.4,
+    # if using 7.0 need to link to: ../extern/lib/${WINDIR}/microsoft/msvc70
+    SET(MATLAB_LIBRARIES_DIR "${MATLAB_ROOT}/extern/lib/${WINDIR}/microsoft/bcc54")
+    # SET(MATLAB_LIBRARIES_DIR "${MATLAB_ROOT}/extern/lib/${WINDIR}/microsoft/bcc50")
+    # SET(MATLAB_LIBRARIES_DIR "${MATLAB_ROOT}/extern/lib/${WINDIR}/microsoft/bcc51")
+  ELSEIF(${CMAKE_GENERATOR} MATCHES "Visual Studio*")
+    # If the compiler is Visual Studio, but not any of the specific
+    # versions above, we try our luck with the microsoft directory
+    SET(MATLAB_LIBRARIES_DIR "${MATLAB_ROOT}/extern/lib/${WINDIR}/microsoft/")
+  ELSE(${CMAKE_GENERATOR} MATCHES "Visual Studio 6")
+    MESSAGE(FATAL_ERROR "Generator not compatible: ${CMAKE_GENERATOR}")
+  ENDIF(${CMAKE_GENERATOR} MATCHES "Visual Studio 6")
+
+  # Get paths to the Matlab MEX libraries
+  FIND_LIBRARY(MATLAB_MEX_LIBRARY
+    libmex
+    ${MATLAB_LIBRARIES_DIR}
+    )
+  FIND_LIBRARY(MATLAB_MX_LIBRARY
+    libmx
+    ${MATLAB_LIBRARIES_DIR}
+    )
+  FIND_LIBRARY(MATLAB_MAT_LIBRARY
+    libmat
+    ${MATLAB_LIBRARIES_DIR}
+    )
+  FIND_LIBRARY(MATLAB_ENG_LIBRARY
+    libeng
+    ${MATLAB_LIBRARIES_DIR}
+    )
+
+  # Get path to the include directory
+  FIND_PATH(MATLAB_INCLUDE_DIR
+    "mex.h"
+    "${MATLAB_ROOT}/extern/include"
+    )
+
+ELSE(WIN32)
+
+  IF((NOT DEFINED MATLAB_ROOT) 
+      OR ("${MATLAB_ROOT}" STREQUAL ""))
+    # get path to the Matlab root directory
+    EXECUTE_PROCESS(
+      COMMAND which matlab
+      COMMAND xargs readlink
+      COMMAND xargs dirname
+      COMMAND xargs dirname
+      COMMAND xargs echo -n
+      OUTPUT_VARIABLE MATLAB_ROOT
+      )
+  ENDIF((NOT DEFINED MATLAB_ROOT) 
+    OR ("${MATLAB_ROOT}" STREQUAL ""))
+
+  # Check if this is a Mac
+  IF(${CMAKE_SYSTEM_NAME} MATCHES "Darwin")
 
-FIND_LIBRARY(MATLAB_LIBRARY2 mx
-  PATHS
-    /Applications/MATLAB_R2012b.app/bin/maci64/
-    /Applications/MATLAB_R2013a.app/bin/maci64/
-    /Applications/MATLAB_R2013b.app/bin/maci64/
-    /Applications/MATLAB_R2014a.app/bin/maci64/
-    /Applications/MATLAB_R2014b.app/bin/maci64/
-  PATH_SUFFIXES
-    a
-    lib64
-    lib
-    dylib
-    NO_DEFAULT_PATH
+    SET(LIBRARY_EXTENSION .dylib)
+
+    # If this is a Mac and the attempts to find MATLAB_ROOT have so far failed, 
+    # we look in the applications folder
+    IF((NOT DEFINED MATLAB_ROOT) OR ("${MATLAB_ROOT}" STREQUAL ""))
+
+    # Search for a version of Matlab available, starting from the most modern one to older versions
+      FOREACH(MATVER "R2014a" "R2014a" "R2013b" "R2013a" "R2012b" "R2012a" "R2011b" "R2011a" "R2010b" "R2010a" "R2009b" "R2009a" "R2008b")
+        IF((NOT DEFINED MATLAB_ROOT) OR ("${MATLAB_ROOT}" STREQUAL ""))
+          IF(EXISTS /Applications/MATLAB_${MATVER}.app)
+            SET(MATLAB_ROOT /Applications/MATLAB_${MATVER}.app)
+    
+          ENDIF(EXISTS /Applications/MATLAB_${MATVER}.app)
+        ENDIF((NOT DEFINED MATLAB_ROOT) OR ("${MATLAB_ROOT}" STREQUAL ""))
+      ENDFOREACH(MATVER)
+
+    ENDIF((NOT DEFINED MATLAB_ROOT) OR ("${MATLAB_ROOT}" STREQUAL ""))
+
+  ELSE(${CMAKE_SYSTEM_NAME} MATCHES "Darwin")
+    SET(LIBRARY_EXTENSION .so)
+
+  ENDIF(${CMAKE_SYSTEM_NAME} MATCHES "Darwin")
+
+  # Get path to the MEX libraries
+  EXECUTE_PROCESS(
+    #COMMAND find "${MATLAB_ROOT}/extern/lib" -name libmex${LIBRARY_EXTENSION} # Peter
+	COMMAND find "${MATLAB_ROOT}/bin" -name libmex${LIBRARY_EXTENSION} # Standard
+    COMMAND xargs echo -n
+    OUTPUT_VARIABLE MATLAB_MEX_LIBRARY
+    )
+  EXECUTE_PROCESS(
+    #COMMAND find "${MATLAB_ROOT}/extern/lib" -name libmx${LIBRARY_EXTENSION} # Peter
+	COMMAND find "${MATLAB_ROOT}/bin" -name libmx${LIBRARY_EXTENSION} # Standard
+    COMMAND xargs echo -n
+    OUTPUT_VARIABLE MATLAB_MX_LIBRARY
+    )
+  EXECUTE_PROCESS(
+    #COMMAND find "${MATLAB_ROOT}/extern/lib" -name libmat${LIBRARY_EXTENSION} # Peter
+	COMMAND find "${MATLAB_ROOT}/bin" -name libmat${LIBRARY_EXTENSION} # Standard
+    COMMAND xargs echo -n
+    OUTPUT_VARIABLE MATLAB_MAT_LIBRARY
+    )
+  EXECUTE_PROCESS(
+    #COMMAND find "${MATLAB_ROOT}/extern/lib" -name libeng${LIBRARY_EXTENSION} # Peter
+	COMMAND find "${MATLAB_ROOT}/bin" -name libeng${LIBRARY_EXTENSION} # Standard
+    COMMAND xargs echo -n
+    OUTPUT_VARIABLE MATLAB_ENG_LIBRARY
+    )
+
+  # Get path to the include directory
+  FIND_PATH(MATLAB_INCLUDE_DIR
+    "mex.h"
+    PATHS "${MATLAB_ROOT}/extern/include"
+    )
+
+ENDIF(WIN32)
+
+# This is common to UNIX and Win32:
+SET(MATLAB_LIBRARIES
+  ${MATLAB_MAT_LIBRARY}
+  ${MATLAB_MEX_LIBRARY}
+  ${MATLAB_MX_LIBRARY}
+  ${MATLAB_ENG_LIBRARY}
 )
 
-if(MATLAB_INCLUDE_DIR AND MATLAB_LIBRARY1 AND MATLAB_LIBRARY2)
-	message(STATUS "Found MATLAB: ${MATLAB_INCLUDE_DIR}")
-  set(MATLAB_LIBRARIES ${MATLAB_LIBRARY1} ${MATLAB_LIBRARY2})
-else(MATLAB_INCLUDE_DIR AND MATLAB_LIBRARY1 AND MATLAB_LIBRARY2)
-	message(FATAL_ERROR "could NOT find MATLAB")
-endif(MATLAB_INCLUDE_DIR AND MATLAB_LIBRARY1 AND MATLAB_LIBRARY2)
+IF(MATLAB_INCLUDE_DIR AND MATLAB_LIBRARIES)
+  SET(MATLAB_FOUND 1)
+ENDIF(MATLAB_INCLUDE_DIR AND MATLAB_LIBRARIES)
+
+MARK_AS_ADVANCED(
+  MATLAB_LIBRARIES
+  MATLAB_MEX_LIBRARY
+  MATLAB_MX_LIBRARY
+  MATLAB_ENG_LIBRARY
+  MATLAB_INCLUDE_DIR
+  MATLAB_FOUND
+  MATLAB_ROOT
+)

tutorial/images/camelhead-laplace-equation.jpg.REMOVED.git-id

@@ -0,0 +1 @@
+39ddde658cbb1e577caeb58083724a81f796956e

tutorial/readme.md

@@ -42,6 +42,8 @@ of these lecture notes links to a cross-platform example application.
     * [301 Slice](#slice)
     * [301 Sort](#sort)
         * [Other Matlab-style functions](#othermatlab-stylefunctions) 
+    * [301 Laplace Equation](#laplaceequation)
+        * [Quadratic energy minimization](#quadraticenergyminimization)
 
 # Compilation Instructions
 
@@ -493,8 +495,140 @@ functionality as common matlab functions.
 - `igl::median` Compute the median per column
 - `igl::mode` Compute the mode per column
 - `igl::orth` Orthogonalization of a basis
+- `igl::setdiff` Set difference of matrix elements
 - `igl::speye` Identity as sparse matrix
 
+## Laplace Equation
+A common linear system in geometry processing is the Laplace equation:
+
+ $∆z = 0$
+
+subject to some boundary conditions, for example Dirichlet boundary conditions
+(fixed value):
+
+ $\left.z\right|_{\partial{S}} = z_{bc}$
+
+In the discrete setting, this begins with the linear system:
+
+ $\mathbf{L} \mathbf{z} = \mathbf{0}$
+
+where $\mathbf{L}$ is the $n \times n$ discrete Laplacian and $\mathbf{z}$ is a
+vector of per-vertex values. Most of $\mathbf{z}$ correspond to interior
+vertices and are unknown, but some of $\mathbf{z}$ represent values at boundary
+vertices. Their values are known so we may move their corresponding terms to
+the right-hand side.
+
+Conceptually, this is very easy if we have sorted $\mathbf{z}$ so that interior
+vertices come first and then boundary vertices:
+
+ $$\left(\begin{array}{cc}
+ \mathbf{L}_{in,in} & \mathbf{L}_{in,b}\\
+ \mathbf{L}_{b,in} & \mathbf{L}_{b,b}\end{array}\right) 
+ \left(\begin{array}{c}
+ \mathbf{z}_{in}\\
+ \mathbf{L}_{b}\end{array}\right) = 
+ \left(\begin{array}{c}
+ \mathbf{0}_{in}\\
+ \mathbf{*}_{b}\end{array}\right)$$ 
+
+The bottom block of equations is no longer meaningful so we'll only consider
+the top block:
+
+ $$\left(\begin{array}{cc}
+ \mathbf{L}_{in,in} & \mathbf{L}_{in,b}\end{array}\right) 
+ \left(\begin{array}{c}
+ \mathbf{z}_{in}\\
+ \mathbf{z}_{b}\end{array}\right) = 
+ \mathbf{0}_{in}$$
+
+Where now we can move known values to the right-hand side:
+
+ $$\mathbf{L}_{in,in} 
+ \mathbf{z}_{in} = -
+ \mathbf{L}_{in,b}
+ \mathbf{z}_{b}$$
+
+Finally we can solve this equation for the unknown values at interior vertices
+$\mathbf{z}_{in}$.
+
+However, probably our vertices are not sorted. One option would be to sort `V`,
+then proceed as above and then _unsort_ the solution `Z` to match `V`. However,
+this solution is not very general.
+
+With array slicing no explicit sort is needed. Instead we can _slice-out_
+submatrix blocks ($\mathbf{L}_{in,in}$, $\mathbf{L}_{in,b}$, etc.) and follow
+the linear algebra above directly. Then we can slice the solution _into_ the
+rows of `Z` corresponding to the interior vertices.
+
+![The `LaplaceEquation` example solves a Laplace equation with Dirichlet
+boundary conditions.](images/camelhead-laplace-equation.jpg)
+
+### Quadratic energy minimization
+
+The same Laplace equation may be equivalently derived by minimizing Dirichlet
+energy subject to the same boundary conditions:
+
+ $\mathop{\text{minimize }}_z \int\limits_S \|\nabla z\|^2 dA$
+
+On our discrete mesh, recall that this becomes
+
+ $\mathop{\text{minimize }}_\mathbf{z}  \mathbf{z}^T \mathbf{G}^T \mathbf{D}
+ \mathbf{G} \mathbf{z} \rightarrow \mathop{\text{minimize }}_\mathbf{z} \mathbf{z}^T \mathbf{L} \mathbf{z}$
+
+The general problem of minimizing some energy over a mesh subject to fixed
+value boundary conditions is so wide spread that libigl has a dedicated api for
+solving such systems. 
+
+Let's consider a general quadratic minimization problem subject to different
+common constraints:
+
+ $$\mathop{\text{minimize }}_\mathbf{z}  \mathbf{z}^T \mathbf{Q} \mathbf{z} +
+ \mathbf{z}^T \mathbf{B} + \text{constant},$$
+
+ subject to
+ 
+ $$\mathbf{z}_b = \mathbf{z}_{bc} \text{ and } \mathbf{A}_{eq} \mathbf{z} =
+ \mathbf{B}_{eq},$$
+
+where 
+
+  - $\mathbf{Q}$ is a (usually sparse) $n \times n$ positive semi-definite
+    matrix of quadratic coefficients (Hessian), 
+  - $\mathbf{B}$ is a $n \times 1$ vector of linear coefficients, 
+  - $\mathbf{z}_b$ is a $|b| \times 1$ portion of
+$\mathbf{z}$ corresponding to boundary or _fixed_ vertices,
+  - $\mathbf{z}_{bc}$ is a $|b| \times 1$ vector of known values corresponding to
+    $\mathbf{z}_b$,
+  - $\mathbf{A}_{eq}$ is a (usually sparse) $m \times n$ matrix of linear
+    equality constraint coefficients (one row per constraint), and
+  - $\mathbf{B}_{eq}$ is a $m \times 1$ vector of linear equality constraint
+    right-hand side values.
+
+This specification is overly general as we could write $\mathbf{z}_b =
+\mathbf{z}_{bc}$ as rows of $\mathbf{A}_{eq} \mathbf{z} =
+\mathbf{B}_{eq}$, but these fixed value constraints appear so often that they
+merit a dedicated place in the API.
+
+In libigl, solving such quadratic optimization problems is split into two
+routines: precomputation and solve. Precomputation only depends on the
+quadratic coefficients, known value indices and linear constraint coefficients:
+
+```cpp
+igl::min_quad_with_fixed_data mqwf;
+igl::min_quad_with_fixed_precompute(Q,b,Aeq,true,mqwf);
+```
+
+The output is a struct `mqwf` which contains the system matrix factorization
+and is used during solving with arbitrary linear terms, known values, and
+constraint right-hand sides:
+
+```cpp
+igl::min_quad_with_fixed_solve(mqwf,B,bc,Beq,Z);
+```
+
+The output `Z` is a $n \times 1$ vector of solutions with fixed values
+correctly placed to match the mesh vertices `V`.
+
 
 [#meyer_2003]: Mark Meyer and Mathieu Desbrun and Peter Schröder and Alan H.  Barr,
 "Discrete Differential-Geometry Operators for Triangulated