support for higher-order biharmonic coords

Former-commit-id: cc7e9e9dcec84f0fe426206aa9a750fa1a62c6a0

include/igl/biharmonic_coordinates.cpp

@@ -23,6 +23,21 @@ IGL_INLINE bool igl::biharmonic_coordinates(
   const Eigen::PlainObjectBase<DerivedT> & T,
   const std::vector<std::vector<SType> > & S,
   Eigen::PlainObjectBase<DerivedW> & W)
+{
+  return biharmonic_coordinates(V,T,S,2,W);
+}
+
+template <
+  typename DerivedV,
+  typename DerivedT,
+  typename SType,
+  typename DerivedW>
+IGL_INLINE bool igl::biharmonic_coordinates(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedT> & T,
+  const std::vector<std::vector<SType> > & S,
+  const int k,
+  Eigen::PlainObjectBase<DerivedW> & W)
 {
   using namespace Eigen;
   using namespace std;
@@ -58,9 +73,23 @@ IGL_INLINE bool igl::biharmonic_coordinates(
     cotmatrix(V,T,L);
     K = N+L;
     massmatrix(V,T,MASSMATRIX_TYPE_DEFAULT,M);
+    // normalize
+    M /= ((VectorXd)M.diagonal()).array().abs().maxCoeff();
     DiagonalMatrix<double,Dynamic> Minv = 
       ((VectorXd)M.diagonal().array().inverse()).asDiagonal();
-    A = K.transpose() * (Minv * K);
+    switch(k)
+    {
+      default:
+        assert(false && "unsupported");
+      case 2:
+        // For C1 smoothness in 2D, one should use bi-harmonic
+        A = K.transpose() * (Minv * K);
+        break;
+      case 3:
+        // For C1 smoothness in 3D, one should use tri-harmonic
+        A = K.transpose() * (Minv * (-L * (Minv * K)));
+        break;
+    }
   }
   // Vertices in point handles
   const size_t mp = 

include/igl/biharmonic_coordinates.h

@@ -70,6 +70,19 @@ namespace igl
     const Eigen::PlainObjectBase<DerivedT> & T,
     const std::vector<std::vector<SType> > & S,
     Eigen::PlainObjectBase<DerivedW> & W);
+  // k  2-->biharmonic, 3-->triharmonic
+  template <
+    typename DerivedV,
+    typename DerivedT,
+    typename SType,
+    typename DerivedW>
+  IGL_INLINE bool biharmonic_coordinates(
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::PlainObjectBase<DerivedT> & T,
+    const std::vector<std::vector<SType> > & S,
+    const int k,
+    Eigen::PlainObjectBase<DerivedW> & W);
+
 };
 #  ifndef IGL_STATIC_LIBRARY
 #    include "biharmonic_coordinates.cpp"