fixed optimization flags in some Makefiles, tetgen wraper, launch medit, mosek wrapper

Former-commit-id: 68bbcbc66e895af55204bb60a38685be93abc43a

Makefile

@@ -6,8 +6,9 @@ include Makefile.conf
 
 # optimized default settings
 all: LFLAGS +=
-all: CFLAGS += -O3 -DNDEBUG -j 
-debug: CFLAGS += -g -Wall -Werror
+OPTFLAGS=-O3 -DNDEBUG -j
+debug: OPTFLAGS= -g -Wall -Werror
+CFLAGS += $(OPTFLAGS)
 
 EXTRA_DIRS=
 ifeq ($(IGL_WITH_TETGEN),1)
@@ -18,6 +19,10 @@ ifeq ($(IGL_WITH_MATLAB),1)
 	# append matlab extra dir
 	EXTRA_DIRS+=include/igl/matlab
 endif
+ifeq ($(IGL_WITH_MOSEK),1)
+	# append mosek extra dir
+	EXTRA_DIRS+=include/igl/mosek
+endif
 
 .PHONY: examples
 .PHONY: extras

Makefile.conf

@@ -1,5 +1,6 @@
 IGL_WITH_TETGEN=1
 IGL_WITH_MATLAB=1
+IGL_WITH_MOSEK=1
 
 #############################################################################
 # FLAGS 

include/igl/faces_first.cpp

@@ -11,6 +11,8 @@ IGL_INLINE void igl::faces_first(
   MatF & RF, 
   VecI & IM)
 {
+  assert(&V != &RV);
+  assert(&F != &RF);
   using namespace std;
   using namespace Eigen;
   vector<bool> in_face(V.rows());

include/igl/faces_first.h

@@ -24,6 +24,13 @@ namespace igl
   //  IM  # faces by 1 list of indices such that: RF = IM(F) and RT = IM(T)
   //    and RV(IM,:) = V
   //
+  //
+  // Example:
+  //   // Tet mesh in (V,T,F)
+  //   faces_first(V,F,IM);
+  //   T = T.unaryExpr(bind1st(mem_fun( static_cast<VectorXi::Scalar&
+  //     (VectorXi::*)(VectorXi::Index)>(&VectorXi::operator())),
+  //     &IM)).eval();
   template <typename MatV, typename MatF, typename VecI>
   IGL_INLINE void faces_first(
     const MatV & V, 

include/igl/launch_medit.cpp

@@ -0,0 +1,56 @@
+#include "launch_medit.h"
+#include "writeMESH.h"
+#include <iostream>
+#include <string>
+#include <sstream>
+
+#define MEDIT_PATH "/opt/local/bin/medit"
+#define TEMP_MESH_FILE "/var/tmp/temp.mesh"
+#define TEMP_MEDIT_FILE "/var/tmp/temp.medit"
+
+template <typename DerivedV, typename DerivedT, typename DerivedF>
+IGL_INLINE int igl::launch_medit(
+  const Eigen::MatrixBase<DerivedV> & V, 
+  const Eigen::MatrixBase<DerivedT> & T,
+  const Eigen::MatrixBase<DerivedF> & F)
+{
+  using namespace std;
+  // Build medit command, end with & so command returns without waiting
+  stringstream command; 
+  command<<MEDIT_PATH<<" "<<TEMP_MESH_FILE<<" "<<TEMP_MEDIT_FILE<<" &"<<endl;
+  bool mesh_saved = writeMESH(TEMP_MESH_FILE,V,T,F);
+  if(!mesh_saved)
+  {
+    return -1;
+  }
+  // Write default medit options
+  const string default_medit_file_contents = 
+    "BackgroundColor 1 1 1\n"
+    "LineColor 0 0 0\n"
+    "WindowSize 1024 800\n"
+    "RenderMode shading + lines\n";
+  FILE * fp = fopen(TEMP_MEDIT_FILE,"w");
+  if(fp == NULL)
+  {
+    cerr<<"^"<<__FUNCTION__<<": Could not write to "<<TEMP_MEDIT_FILE<<endl;
+    return -1;
+  }
+  fprintf(fp,"%s",default_medit_file_contents.c_str());
+  fclose(fp);
+
+  try
+  {
+    return system(command.str().c_str());
+  }catch(int e)
+  {
+    cerr<<"^"<<__FUNCTION__<<": Calling to medit crashed..."<<endl;
+    return -1;
+  }
+  // Could clean up and delete these files but not really worth it
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+template int igl::launch_medit<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&);
+#endif
+

include/igl/launch_medit.h

@@ -0,0 +1,35 @@
+#ifndef IGL_LAUNCH_MEDIT_H
+#define IGL_LAUNCH_MEDIT_H
+#include "igl_inline.h"
+
+#include <Eigen/Dense>
+
+namespace igl 
+{
+  // Writes the tetmesh in (V,T,F) to a temporary file, opens it with medit
+  // (forking with a system call) and returns
+  //
+  //
+  // Templates:
+  //   DerivedV  real-value: i.e. from MatrixXd
+  //   DerivedT  integer-value: i.e. from MatrixXi
+  //   DerivedF  integer-value: i.e. from MatrixXi
+  // Inputs:
+  //   V  double matrix of vertex positions  #V by 3
+  //   T  #T list of tet indices into vertex positions
+  //   F  #F list of face indices into vertex positions
+  // Returns returned value of system call (probably not useful because of the
+  //   fork)
+  template <typename DerivedV, typename DerivedT, typename DerivedF>
+  IGL_INLINE int launch_medit(
+    const Eigen::MatrixBase<DerivedV> & V, 
+    const Eigen::MatrixBase<DerivedT> & T,
+    const Eigen::MatrixBase<DerivedF> & F);
+}
+
+#ifdef IGL_HEADER_ONLY
+#  include "launch_medit.cpp"
+#endif
+
+#endif
+

include/igl/matlab/Makefile

@@ -1,10 +1,12 @@
 include ../../../Makefile.conf
+all: CFLAGS += -O3 -DNDEBUG -j 
+debug: CFLAGS += -g -Wall -Werror
 
 .PHONY: all
-all: libmatlab
+all: libiglmatlab
 
 .PHONY: libmatlab
-libmatlab: obj ../../../lib/libiglmatlab.a
+libiglmatlab: obj ../../../lib/libiglmatlab.a
 
 CPP_FILES=$(wildcard *.cpp)
 OBJ_FILES=$(addprefix obj/,$(notdir $(CPP_FILES:.cpp=.o)))

include/igl/mosek/Makefile

@@ -0,0 +1,41 @@
+include ../../../Makefile.conf
+
+.PHONY: all
+all: libiglmosek
+
+.PHONY: libiglmosek
+libiglmosek: obj ../../../lib/libiglmosek.a
+
+CPP_FILES=$(wildcard *.cpp)
+OBJ_FILES=$(addprefix obj/,$(notdir $(CPP_FILES:.cpp=.o)))
+
+# include igl headers
+INC+=-I../../../include/
+
+# EXPECTS THAT CFLAGS IS ALREADY SET APPROPRIATELY 
+
+# Eigen dependency
+EIGEN3_INC=-I/opt/local/include/eigen3 -I/opt/local/include/eigen3/unsupported
+INC+=$(EIGEN3_INC)
+
+# mosek dependency
+# TODO: linux, 32 bit etc
+MOSEKPLATFORM=osx64x86
+MOSEK=/usr/local/mosek
+MOSEK_INC=-I$(MOSEK)/6/tools/platform/$(MOSEKPLATFORM)/h
+MOSEK_LIB=-L$(MOSEK)/6/tools/platform/$(MOSEKPLATFORM)/bin -lmosek64
+INC+=$(MOSEK_INC)
+
+obj: 
+	mkdir -p obj
+
+../../../lib/libiglmosek.a: $(OBJ_FILES)
+	rm -f $@
+	ar cqs $@ $(OBJ_FILES)
+
+obj/%.o: %.cpp %.h
+	g++ $(AFLAGS) $(CFLAGS) -c -o $@ $< $(INC)
+
+clean:
+	rm -f obj/*.o
+	rm -f ../../../lib/libiglmosek.a

include/igl/mosek/mosek_guarded.cpp

@@ -0,0 +1,17 @@
+#include "mosek_guarded.h"
+#include <iostream>
+
+IGL_INLINE MSKrescodee igl::mosek_guarded(const MSKrescodee r)
+{
+  using namespace std;
+  if(r != MSK_RES_OK)
+  {
+    /* In case of an error print error code and description. */      
+    char symname[MSK_MAX_STR_LEN];
+    char desc[MSK_MAX_STR_LEN];
+    MSK_getcodedesc(r,symname,desc);
+    cerr<<"MOSEK ERROR ("<<r<<"): "<<symname<<" - '"<<desc<<"'"<<endl;
+  }
+  return r;
+}
+

include/igl/mosek/mosek_guarded.h

@@ -0,0 +1,21 @@
+#ifndef IGL_MOSEK_GUARDED_H
+#define IGL_MOSEK_GUARDED_H
+#include "../igl_inline.h"
+
+#include "mosek.h"
+namespace igl
+{
+  // Little function to wrap around mosek call to handle errors
+  // 
+  // Inputs:
+  //   r  mosek error code returned from mosek call
+  // Returns r untouched
+  IGL_INLINE MSKrescodee mosek_guarded(const MSKrescodee r);
+}
+
+#ifdef IGL_HEADER_ONLY
+#  include "mosek_guarded.cpp"
+#endif
+
+#endif
+

include/igl/mosek/mosek_quadprog.cpp

@@ -0,0 +1,226 @@
+#include "mosek_quadprog.h"
+#include "mosek_guarded.h"
+#include <cstdio>
+#include "../verbose.h"
+
+// Little function mosek needs in order to know how to print to std out
+static void MSKAPI printstr(void *handle, char str[])
+{
+  printf("%s",str);
+}
+
+template <typename Index, typename Scalar>
+IGL_INLINE bool igl::mosek_quadprog(
+  const Index n,
+  const std::vector<Index> & Qi,
+  const std::vector<Index> & Qj,
+  const std::vector<Scalar> & Qv,
+  const std::vector<Scalar> & c,
+  const Scalar cf,
+  const Index m,
+  const std::vector<Index> & Ari,
+  const std::vector<Index> & Acp,
+  const std::vector<Scalar> & Av,
+  const std::vector<Scalar> & lc,
+  const std::vector<Scalar> & uc,
+  const std::vector<Scalar> & lx,
+  const std::vector<Scalar> & ux,
+  std::vector<Scalar> & x)
+{
+  // I J V vectors of Q should all be same length
+  assert(Qv.size() == Qi.size());
+  assert(Qv.size() == Qj.size());
+  // number of columns in linear constraint matrix must be ≤ number of
+  // variables
+  assert(Acp.size() == (n+1));
+  // linear bound vectors must be size of number of constraints or empty
+  assert( (lc.size() == m) || (lc.size() == 0));
+  assert( (uc.size() == m) || (uc.size() == 0));
+  // constant bound vectors must be size of number of variables or empty
+  assert( (lx.size() == n) || (lx.size() == 0));
+  assert( (ux.size() == n) || (ux.size() == 0));
+
+  // allocate space for solution in x
+  x.resize(n);
+
+  // variables for mosek task, env and result code
+  MSKenv_t env;
+  MSKtask_t task;
+
+  // Create the MOSEK environment
+  mosek_guarded(MSK_makeenv(&env,NULL,NULL,NULL,NULL));
+  /* Directs the log stream to the 'printstr' function. */
+  mosek_guarded(MSK_linkfunctoenvstream(env,MSK_STREAM_LOG,NULL,printstr));
+  // initialize mosek environment
+  mosek_guarded(MSK_initenv(env));
+  // Create the optimization task
+  mosek_guarded(MSK_maketask(env,m,n,&task));
+  verbose("Creating task with %ld linear constraints and %ld variables...\n",m,n);
+  // Tell mosek how to print to std out
+  mosek_guarded(MSK_linkfunctotaskstream(task,MSK_STREAM_LOG,NULL,printstr));
+  // Give estimate of number of variables
+  mosek_guarded(MSK_putmaxnumvar(task,n));
+  if(m>0)
+  {
+    // Give estimate of number of constraints
+    mosek_guarded(MSK_putmaxnumcon(task,m));
+    // Give estimate of number of non zeros in A
+    mosek_guarded(MSK_putmaxnumanz(task,Av.size()));
+  }
+  // Give estimate of number of non zeros in Q
+  mosek_guarded(MSK_putmaxnumqnz(task,Qv.size()));
+  if(m>0)
+  {
+    // Append 'm' empty constraints, the constrainst will initially have no
+    // bounds
+    mosek_guarded(MSK_append(task,MSK_ACC_CON,m));
+  }
+  // Append 'n' variables
+  mosek_guarded(MSK_append(task,MSK_ACC_VAR,n));
+  // add a contant term to the objective
+  mosek_guarded(MSK_putcfix(task,cf));
+
+  // loop over variables
+  for(int j = 0;j<n;j++)
+  {
+    if(c.size() > 0)
+    {
+      // Set linear term c_j in the objective
+      mosek_guarded(MSK_putcj(task,j,c[j]));
+    }
+
+    // Set constant bounds on variable j
+    if(lx[j] == ux[j])
+    {
+      mosek_guarded(MSK_putbound(task,MSK_ACC_VAR,j,MSK_BK_FX,lx[j],ux[j]));
+    }else
+    {
+      mosek_guarded(MSK_putbound(task,MSK_ACC_VAR,j,MSK_BK_RA,lx[j],ux[j]));
+    }
+
+    if(m>0)
+    {
+      // Input column j of A
+      mosek_guarded(
+        MSK_putavec(
+          task,
+          MSK_ACC_VAR,
+          j,
+          Acp[j+1]-Acp[j],
+          &Ari[Acp[j]],
+          &Av[Acp[j]])
+        );
+    }
+  }
+
+  // loop over constraints
+  for(size_t i = 0;i<m;i++)
+  {
+    // put bounds  on constraints
+    mosek_guarded(MSK_putbound(task,MSK_ACC_CON,i,MSK_BK_RA,lc[i],uc[i]));
+  }
+
+  // Input Q for the objective (REMEMBER Q SHOULD ONLY BE LOWER TRIANGLE)
+  mosek_guarded(MSK_putqobj(task,Qv.size(),&Qi[0],&Qj[0],&Qv[0]));
+
+  //// Set up task parameters
+
+  // set tolerance
+  //mosek_guarded(
+  //  MSK_putdouparam(task,MSK_DPAR_INTPNT_TOL_DFEAS,1e-8));
+  //mosek_guarded(
+  //  MSK_putdouparam(task,MSK_DPAR_INTPNT_TOL_DSAFE,1.0));
+  // 1.0 means optimizer is very leniant about declaring model infeasible
+  //mosek_guarded(
+  //  MSK_putdouparam(task,MSK_DPAR_INTPNT_TOL_INFEAS,1e-8));
+  //mosek_guarded(
+  //  MSK_putdouparam(task,MSK_DPAR_INTPNT_TOL_PATH,1e-8));
+  //mosek_guarded(
+  //  MSK_putdouparam(task,MSK_DPAR_INTPNT_TOL_PFEAS,1e-8));
+
+  // Hard to say if this is doing anything, probably nothing dramatic
+  mosek_guarded(MSK_putdouparam(task,MSK_DPAR_INTPNT_TOL_PSAFE,1e2));
+
+  // >1e0 NONSOLUTION
+  // 1e-1 artifacts in deformation
+  // 1e-3 artifacts in isolines
+  // 1e-4 seems safe
+  // 1e-8 MOSEK DEFAULT SOLUTION
+  mosek_guarded(MSK_putdouparam(task,MSK_DPAR_INTPNT_TOL_REL_GAP,1e-8));
+  //mosek_guarded(MSK_putdouparam(task,MSK_DPAR_INTPNT_TOL_REL_STEP,0.9999));
+
+  // Turn off presolving
+  //mosek_guarded(
+  //  MSK_putintparam(task,MSK_IPAR_PRESOLVE_USE,MSK_PRESOLVE_MODE_OFF));
+
+  // Force particular matrix reordering method
+  // MSK_ORDER_METHOD_NONE cuts time in half roughly, since half the time is
+  //   usually spent reordering the matrix
+  // !! WARNING Setting this parameter to anything but MSK_ORDER_METHOD_FREE
+  //   seems to have the effect of setting it to MSK_ORDER_METHOD_NONE
+  //   *Or maybe Mosek is spending a bunch of time analyzing the matrix to
+  //   choose the right ordering method when really any of them are
+  //   instantaneous
+  mosek_guarded(
+    MSK_putintparam(task,MSK_IPAR_INTPNT_ORDER_METHOD,MSK_ORDER_METHOD_NONE));
+
+  // Turn off convexity check
+  mosek_guarded(
+    MSK_putintparam(task,MSK_IPAR_CHECK_CONVEXITY,MSK_CHECK_CONVEXITY_NONE));
+
+  // Force using multiple threads, not sure if MOSEK is properly destorying
+  //extra threads...
+  mosek_guarded(MSK_putintparam(task,MSK_IPAR_INTPNT_NUM_THREADS,6));
+  
+  // Force turn off data check
+  mosek_guarded(MSK_putintparam(task,MSK_IPAR_DATA_CHECK,MSK_OFF));
+
+  // Now the optimizer has been prepared
+  MSKrescodee trmcode;
+  // run the optimizer
+  mosek_guarded(MSK_optimizetrm(task,&trmcode));
+
+  // Print a summary containing information about the solution for debugging
+  // purposes
+  MSK_solutionsummary(task,MSK_STREAM_LOG);
+
+  // Get status of solution
+  MSKsolstae solsta;
+  MSK_getsolutionstatus(task,MSK_SOL_ITR,NULL,&solsta);
+
+  bool success = false;
+  switch(solsta)
+  {
+    case MSK_SOL_STA_OPTIMAL:   
+    case MSK_SOL_STA_NEAR_OPTIMAL:
+      MSK_getsolutionslice(task,MSK_SOL_ITR,MSK_SOL_ITEM_XX,0,n,&x[0]);
+      printf("Optimal primal solution\n");
+      //for(size_t j=0; j<n; ++j)
+      //{
+      //  printf("x[%ld]: %g\n",j,x[j]);
+      //}
+      success = true;
+      break;
+    case MSK_SOL_STA_DUAL_INFEAS_CER:
+    case MSK_SOL_STA_PRIM_INFEAS_CER:
+    case MSK_SOL_STA_NEAR_DUAL_INFEAS_CER:
+    case MSK_SOL_STA_NEAR_PRIM_INFEAS_CER:  
+      printf("Primal or dual infeasibility certificate found.\n");
+      break;
+    case MSK_SOL_STA_UNKNOWN:
+      printf("The status of the solution could not be determined.\n");
+      break;
+    default:
+      printf("Other solution status.");
+      break;
+  }
+
+  MSK_deletetask(&task);
+  MSK_deleteenv(&env);
+
+  return success;
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+#endif

include/igl/mosek/mosek_quadprog.h

@@ -0,0 +1,86 @@
+#ifndef IGL_MOSEK_QUADPROG_H
+#define IGL_MOSEK_QUADPROG_H
+#include "../igl_inline.h"
+#include <vector>
+namespace igl
+{
+  // Solve a convex quadratic optimization problem with linear and constant
+  // bounds, that is:
+  //
+  // Minimize: ½ * xT * Q⁰ * x + cT * x + cf
+  //
+  // Subject to: lc ≤ Ax ≤ uc
+  //             lx ≤ x ≤ ux
+  //
+  // where we are trying to find the optimal vector of values x. 
+  //
+  // Note: Q⁰ must be symmetric and the ½ is a convention of MOSEK
+  //
+  // Note: Because of how MOSEK accepts different parts of the system, Q should
+  // be stored in IJV (aka Coordinate) format and should only include entries in
+  // the lower triangle. A should be stored in Column compressed (aka Harwell
+  // Boeing) format. As described:
+  // http://netlib.org/linalg/html_templates/node92.html
+  // or
+  // http://en.wikipedia.org/wiki/Sparse_matrix
+  //   #Compressed_sparse_column_.28CSC_or_CCS.29
+  // 
+  //
+  // Templates:
+  //   Index  type for index variables
+  //   Scalar  type for floating point variables (gets cast to double?)
+  // Input:
+  //   n  number of variables, i.e. size of x
+  //   Qi  vector of qnnz row indices of non-zeros in LOWER TRIANGLE ONLY of Q⁰
+  //   Qj  vector of qnnz column indices of non-zeros in LOWER TRIANGLE ONLY of 
+  //     Q⁰
+  //   Qv  vector of qnnz values of non-zeros in LOWER TRIANGLE ONLY of Q⁰, 
+  //     such that:
+  //     Q⁰(Qi[k],Qj[k]) = Qv[k] for k ∈ [0,Qnnz-1], where Qnnz is the number of
+  //     non-zeros in Q⁰
+  //   c   (optional) vector of n values of c, transpose of coefficient row vector
+  //     of linear terms, EMPTY means c == 0
+  //   cf  (optional) value of constant term in objective, 0 means cf == 0, so
+  //     optional only in the sense that it is mandatory
+  //   m  number of constraints, therefore also number of rows in linear
+  //     constraint coefficient matrix A, and in linear constraint bound vectors 
+  //     lc and uc
+  //   Ari  vector of row indices corresponding to non-zero values of A,
+  //   Acp  vector of indices into Ari and Av of the first entry for each column
+  //     of A, size(Acp) = (# columns of A) + 1 = n + 1
+  //   Av  vector of non-zero values of A, in column compressed order
+  //   lc  vector of m linear constraint lower bounds
+  //   uc  vector of m linear constraint upper bounds
+  //   lx  vector of n constant lower bounds
+  //   ux  vector of n constant upper bounds
+  // Output:
+  //   x  vector of size n to hold output of optimization
+  // Return:
+  //   true only if optimization was successful with no errors
+  // 
+  // Note: All indices are 0-based
+  //
+  template <typename Index, typename Scalar>
+  IGL_INLINE bool mosek_quadprog(
+    const Index n,
+    const std::vector<Index> & Qi,
+    const std::vector<Index> & Qj,
+    const std::vector<Scalar> & Qv,
+    const std::vector<Scalar> & c,
+    const Scalar cf,
+    const Index m,
+    const std::vector<Index> & Ari,
+    const std::vector<Index> & Acp,
+    const std::vector<Scalar> & Av,
+    const std::vector<Scalar> & lc,
+    const std::vector<Scalar> & uc,
+    const std::vector<Scalar> & lx,
+    const std::vector<Scalar> & ux,
+    std::vector<Scalar> & x);
+}
+
+#ifdef IGL_HEADER_ONLY
+#  include "mosek_quadprog.cpp"
+#endif
+
+#endif

include/igl/tetgen/Makefile

@@ -1,10 +1,14 @@
-include ../../../Makefile.conf
 
 .PHONY: all
-all: libtetgen
+all: libigltetgen
+debug: libigltetgen
+
+include ../../../Makefile.conf
+all: CFLAGS += -O3 -DNDEBUG -j 
+debug: CFLAGS += -g -Wall -Werror
 
 .PHONY: libtetgen
-libtetgen: obj ../../../lib/libigltetgen.a
+libigltetgen: obj ../../../lib/libigltetgen.a
 
 CPP_FILES=$(wildcard *.cpp)
 OBJ_FILES=$(addprefix obj/,$(notdir $(CPP_FILES:.cpp=.o)))

include/igl/tetgen/mesh_to_tetgenio.h

@@ -2,6 +2,7 @@
 #define IGL_MESH_TO_TETGENIO_H
 #include "../igl_inline.h"
 
+#define TETLIBRARY 
 #include "tetgen.h" // Defined tetgenio, REAL
 #include <vector>
 #include <Eigen/Core>

include/igl/tetgen/read_into_tetgenio.h

@@ -3,6 +3,7 @@
 #include "../igl_inline.h"
 
 #include <string>
+#define TETLIBRARY 
 #include "tetgen.h" // Defined tetgenio, REAL
 
 namespace igl

include/igl/tetgen/tetgenio_to_tetmesh.cpp

@@ -53,7 +53,7 @@ IGL_INLINE bool igl::tetgenio_to_tetmesh(
   }
   assert(min_index >= 0);
   assert(max_index >= 0);
-  assert(max_index < V.size());
+  assert(max_index < (int)V.size());
 
   return true;
 }

include/igl/tetgen/tetgenio_to_tetmesh.h

@@ -2,6 +2,7 @@
 #define IGL_TETGENIO_TO_TETMESH_H
 #include "../igl_inline.h"
 
+#define TETLIBRARY 
 #include "tetgen.h" // Defined tetgenio, REAL
 #include <vector>
 #include <Eigen/Core>

include/igl/tetgen/tetrahedralize.cpp

@@ -0,0 +1,100 @@
+#include "tetrahedralize.h"
+#include "mesh_to_tetgenio.h"
+#include "tetgenio_to_tetmesh.h"
+
+// IGL includes 
+#include <igl/matrix_to_list.h>
+#include <igl/list_to_matrix.h>
+#include <igl/boundary_faces.h>
+
+// STL includes
+#include <cassert>
+#include <iostream>
+
+IGL_INLINE int igl::tetrahedralize(
+  const std::vector<std::vector<REAL > > & V, 
+  const std::vector<std::vector<int> > & F, 
+  const std::string switches,
+  std::vector<std::vector<REAL > > & TV, 
+  std::vector<std::vector<int > > & TT, 
+  std::vector<std::vector<int> > & TF)
+{
+  using namespace std;
+  tetgenio in,out;
+  bool success;
+  success = mesh_to_tetgenio(V,F,in);
+  if(!success)
+  {
+    return -1;
+  }
+  try
+  {
+    char * cswitches = new char[switches.size() + 1];
+    std::copy(switches.begin(), switches.end(), cswitches);
+    ::tetrahedralize(cswitches,&in, &out);
+    delete[] cswitches;
+  }catch(int e)
+  {
+    cerr<<"^"<<__FUNCTION__<<": TETGEN CRASHED... KABOOOM!!!"<<endl;
+    return 1;
+  }
+  if(out.numberoftetrahedra == 0)
+  {
+    cerr<<"^"<<__FUNCTION__<<": Tetgen failed to create tets"<<endl;
+    return 2;
+  }
+  success = tetgenio_to_tetmesh(out,TV,TT);
+  if(!success)
+  {
+    return -1;
+  }
+  boundary_faces(TT,TF);
+  return 0;
+}
+
+template <
+  typename DerivedV, 
+  typename DerivedF, 
+  typename DerivedTV, 
+  typename DerivedTT, 
+  typename DerivedTF>
+IGL_INLINE int igl::tetrahedralize(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  const std::string switches,
+  Eigen::PlainObjectBase<DerivedTV>& TV,
+  Eigen::PlainObjectBase<DerivedTT>& TT,
+  Eigen::PlainObjectBase<DerivedTF>& TF)
+{
+  using namespace igl;
+  using namespace std;
+  vector<vector<REAL> > vV,vTV;
+  vector<vector<int> > vF,vTT,vTF;
+  matrix_to_list(V,vV);
+  matrix_to_list(F,vF);
+  int e = tetrahedralize(vV,vF,switches,vTV,vTT,vTF);
+  if(e == 0)
+  {
+    bool TV_rect = list_to_matrix(vTV,TV);
+    if(!TV_rect)
+    {
+      return false;
+    }
+    bool TT_rect = list_to_matrix(vTT,TT);
+    if(!TT_rect)
+    {
+      return false;
+    }
+    bool TF_rect = list_to_matrix(vTF,TF);
+    if(!TF_rect)
+    {
+      return false;
+    }
+  }
+  return e;
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+template int igl::tetrahedralize<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >, 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> >&);
+#endif

include/igl/tetgen/tetrahedralize.h

@@ -0,0 +1,65 @@
+#ifndef IGL_TETRAHEDRALIZE_H
+#define IGL_TETRAHEDRALIZE_H
+#include "../igl_inline.h"
+
+#include <vector>
+#include <string>
+#include <Eigen/Core>
+#define TETLIBRARY 
+#include "tetgen.h" // Defined REAL
+
+namespace igl
+{
+  // Mesh the interior of a surface mesh (V,F) using tetgen
+  //
+  // Inputs:
+  //   V  #V by 3 vertex position list
+  //   F  #F list of polygon face indices into V (0-indexed)
+  //   switches  string of tetgen options (See tetgen documentation) e.g.
+  //     "pq1.414a0.01" tries to mesh the interior of a given surface with
+  //       quality and area constraints
+  //     "" will mesh the convex hull constrained to pass through V (ignores F)
+  // Outputs:
+  //   TV  #V by 3 vertex position list
+  //   TT  #T by 4 list of tet face indices
+  //   TF  #F by 3 list of trianlge face indices
+  // Returns status:
+  //   0 success
+  //   1 tetgen threw exception
+  //   2 tetgen did not crash but could not create any tets (probably there are
+  //     holes, duplicate faces etc.)
+  //   -1 other error
+  IGL_INLINE int tetrahedralize(
+    const std::vector<std::vector<REAL > > & V, 
+    const std::vector<std::vector<int> > & F, 
+    const std::string switches,
+    std::vector<std::vector<REAL > > & TV, 
+    std::vector<std::vector<int > > & TT, 
+    std::vector<std::vector<int> > & TF);
+  
+  // Wrapper with Eigen types
+  // Templates:
+  //   DerivedV  real-value: i.e. from MatrixXd
+  //   DerivedF  integer-value: i.e. from MatrixXi
+  template <
+    typename DerivedV, 
+    typename DerivedF, 
+    typename DerivedTV, 
+    typename DerivedTT, 
+    typename DerivedTF>
+  IGL_INLINE int tetrahedralize(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    const std::string switches,
+    Eigen::PlainObjectBase<DerivedTV>& TV,
+    Eigen::PlainObjectBase<DerivedTT>& TT,
+    Eigen::PlainObjectBase<DerivedTF>& TF);
+}
+
+
+#ifdef IGL_HEADER_ONLY
+#  include "tetrahedralize.cpp"
+#endif
+
+#endif
+

include/igl/upsample.cpp

@@ -1,7 +1,6 @@
 #include "upsample.h"
 
 #include "tt.h"
-#include "adjacency_list.h"
 #include <Eigen/Dense>
 
 template <typename MatV, typename MatF>
@@ -40,8 +39,9 @@ IGL_INLINE void igl::upsample( const MatV & V, const MatF & F, MatV & NV, MatF &
   int n_odd = V.rows();
   int n_even = counter;
 
-  Eigen::DynamicSparseMatrix<double> SUBD(V.rows()+n_even,V.rows());
-  SUBD.reserve(15 * (V.rows()+n_even));
+  // Not sure what this is
+  //Eigen::DynamicSparseMatrix<double> SUBD(V.rows()+n_even,V.rows());
+  //SUBD.reserve(15 * (V.rows()+n_even));
   
   // Preallocate NV and NF
   NV = MatV(V.rows()+n_even,V.cols());

include/igl/writeMESH.cpp

@@ -2,26 +2,44 @@
 
 #include <cstdio>
 #include "verbose.h"
+#include "matrix_to_list.h"
+#include <Eigen/Core>
 
 template <typename Scalar, typename Index>
 IGL_INLINE bool igl::writeMESH(
   const std::string mesh_file_name,
-  std::vector<std::vector<Scalar > > & V,
-  std::vector<std::vector<Index > > & T,
-  std::vector<std::vector<Index > > & F)
+  const std::vector<std::vector<Scalar > > & V,
+  const std::vector<std::vector<Index > > & T,
+  const std::vector<std::vector<Index > > & F)
 {
-  // not implemented but should be
-  assert(false);
-  return false;
+  Eigen::MatrixXd mV;
+  Eigen::MatrixXi mT,mF;
+  bool is_rect;
+  is_rect = list_to_matrix(V,mV);
+  if(!is_rect)
+  {
+    return false;
+  }
+  is_rect = list_to_matrix(T,mT);
+  if(!is_rect)
+  {
+    return false;
+  }
+  is_rect = list_to_matrix(F,mF);
+  if(!is_rect)
+  {
+    return false;
+  }
+  return igl::writeMESH(mesh_file_name,mV,mT,mF);
 }
 
-#include <Eigen/Core>
 
+template <typename DerivedV, typename DerivedT, typename DerivedF>
 IGL_INLINE bool igl::writeMESH(
   const std::string str,
-  Eigen::MatrixXd& V,
-  Eigen::MatrixXi& T,
-  Eigen::MatrixXi& F)
+  const Eigen::MatrixBase<DerivedV> & V, 
+  const Eigen::MatrixBase<DerivedT> & T,
+  const Eigen::MatrixBase<DerivedF> & F)
 {
   using namespace std;
   using namespace igl;
@@ -87,4 +105,5 @@ IGL_INLINE bool igl::writeMESH(
 
 #ifndef IGL_HEADER_ONLY
 // Explicit template specialization
+template bool igl::writeMESH<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&);
 #endif

include/igl/writeMESH.h

@@ -15,28 +15,31 @@ namespace igl
   //   Index  type for indices (will be cast to int)
   // Input:
   //   mesh_file_name  path of .mesh file
-  // Outputs:
   //   V  double matrix of vertex positions  #V by 3
   //   T  #T list of tet indices into vertex positions
   //   F  #F list of face indices into vertex positions
   template <typename Scalar, typename Index>
   IGL_INLINE bool writeMESH(
     const std::string mesh_file_name,
-    std::vector<std::vector<Scalar > > & V,
-    std::vector<std::vector<Index > > & T,
-    std::vector<std::vector<Index > > & F);
+    const std::vector<std::vector<Scalar > > & V,
+    const std::vector<std::vector<Index > > & T,
+    const std::vector<std::vector<Index > > & F);
 
+  // Templates:
+  //   DerivedV  real-value: i.e. from MatrixXd
+  //   DerivedT  integer-value: i.e. from MatrixXi
+  //   DerivedF  integer-value: i.e. from MatrixXi
   // Input:
   //   mesh_file_name  path of .mesh file
-  // Outputs:
   //   V  eigen double matrix #V by 3
   //   T  eigen int matrix #T by 4
   //   F  eigen int matrix #F by 3
+  template <typename DerivedV, typename DerivedT, typename DerivedF>
   IGL_INLINE bool writeMESH(
     const std::string str,
-    Eigen::MatrixXd& V,
-    Eigen::MatrixXi& T,
-    Eigen::MatrixXi& F);
+    const Eigen::MatrixBase<DerivedV> & V, 
+    const Eigen::MatrixBase<DerivedT> & T,
+    const Eigen::MatrixBase<DerivedF> & F);
 }
 
 #ifdef IGL_HEADER_ONLY