* added gradMat function

Former-commit-id: 07adecd90cb5ade1f02b829a9ed0eb12e1b27b28

include/igl/gradMat.cpp

@@ -0,0 +1,97 @@
+#include "gradMat.h"
+
+template <typename T, typename S>
+IGL_INLINE void igl::gradMat(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &V,
+  const Eigen::Matrix<S, Eigen::Dynamic, Eigen::Dynamic> &F,
+  Eigen::SparseMatrix<T> &G )
+{
+  Eigen::PlainObjectBase<Eigen::Matrix<T,Eigen::Dynamic,3> > eperp21, eperp13;
+  eperp21.resize(F.rows(),3);
+  eperp13.resize(F.rows(),3);
+
+  for (int i=0;i<F.rows();++i)
+  {
+    // renaming indices of vertices of triangles for convenience
+    int i1 = F(i,0);
+    int i2 = F(i,1);
+    int i3 = F(i,2);
+    
+    // #F x 3 matrices of triangle edge vectors, named after opposite vertices
+    Eigen::Matrix<T, 1, 3> v32 = V.row(i3) - V.row(i2);
+    Eigen::Matrix<T, 1, 3> v13 = V.row(i1) - V.row(i3);
+    Eigen::Matrix<T, 1, 3> v21 = V.row(i2) - V.row(i1);
+    
+    // area of parallelogram is twice area of triangle
+    // area of parallelogram is || v1 x v2 || 
+    Eigen::Matrix<T, 1, 3> n  = v32.cross(v13); 
+    
+    // This does correct l2 norm of rows, so that it contains #F list of twice
+    // triangle areas
+    double dblA = std::sqrt(n.dot(n));
+    
+    // now normalize normals to get unit normals
+    Eigen::Matrix<T, 1, 3> u = n / dblA;
+    
+    // rotate each vector 90 degrees around normal
+    double norm21 = std::sqrt(v21.dot(v21));
+    double norm13 = std::sqrt(v13.dot(v13));
+    eperp21.row(i) = u.cross(v21);
+    eperp21.row(i) = eperp21.row(i) / std::sqrt(eperp21.row(i).dot(eperp21.row(i)));
+    eperp21.row(i) *= norm21 / dblA;
+    eperp13.row(i) = u.cross(v13);
+    eperp13.row(i) = eperp13.row(i) / std::sqrt(eperp13.row(i).dot(eperp13.row(i)));
+    eperp13.row(i) *= norm13 / dblA;
+  }
+
+  std::vector<int> rs;
+  rs.reserve(F.rows()*4*3);
+  std::vector<int> cs;
+  cs.reserve(F.rows()*4*3);
+  std::vector<double> vs;
+  vs.reserve(F.rows()*4*3);
+
+  // row indices
+  for(int r=0;r<3;r++)
+  {
+    for(int j=0;j<4;j++)
+    {
+      for(int i=r*F.rows();i<(r+1)*F.rows();i++) rs.push_back(i);
+    }
+  }
+
+  // column indices
+  for(int r=0;r<3;r++)
+  {
+    for(int i=0;i<F.rows();i++) cs.push_back(F(i,1));
+    for(int i=0;i<F.rows();i++) cs.push_back(F(i,0));
+    for(int i=0;i<F.rows();i++) cs.push_back(F(i,2));
+    for(int i=0;i<F.rows();i++) cs.push_back(F(i,0));
+  }
+  
+  // values
+  for(int i=0;i<F.rows();i++) vs.push_back(eperp13(i,0));
+  for(int i=0;i<F.rows();i++) vs.push_back(-eperp13(i,0));
+  for(int i=0;i<F.rows();i++) vs.push_back(eperp21(i,0));
+  for(int i=0;i<F.rows();i++) vs.push_back(-eperp21(i,0));
+  for(int i=0;i<F.rows();i++) vs.push_back(eperp13(i,1));
+  for(int i=0;i<F.rows();i++) vs.push_back(-eperp13(i,1));
+  for(int i=0;i<F.rows();i++) vs.push_back(eperp21(i,1));
+  for(int i=0;i<F.rows();i++) vs.push_back(-eperp21(i,1));
+  for(int i=0;i<F.rows();i++) vs.push_back(eperp13(i,2));
+  for(int i=0;i<F.rows();i++) vs.push_back(-eperp13(i,2));
+  for(int i=0;i<F.rows();i++) vs.push_back(eperp21(i,2));
+  for(int i=0;i<F.rows();i++) vs.push_back(-eperp21(i,2));
+
+  // create sparse gradient operator matrix
+  G.resize(3*F.rows(),V.rows()); 
+  std::vector<Eigen::Triplet<T> > triplets;
+  for (int i=0;i<vs.size();++i)
+  {
+    triplets.push_back(Eigen::Triplet<T>(rs[i],cs[i],vs[i]));
+  }
+  G.setFromTriplets(triplets.begin(), triplets.end());
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+#endif

include/igl/gradMat.h

@@ -0,0 +1,44 @@
+//
+//  gradMat.h
+//
+//  Created by Christian Schüller on 07/03/13.
+//  Copyright (c) 2013 ETHZ IGL. All rights reserved.
+//
+
+#ifndef IGL_GRAD_MAT_H
+#define IGL_GRAD_MAT_H
+#include "igl_inline.h"
+
+#include <Eigen/Core>
+
+namespace igl {
+  // GRAD
+  // G = grad(V,F)
+  //
+  // Compute the numerical gradient operator
+  //
+  // Inputs:
+  //   V  #vertices by 3 list of mesh vertex positions
+  //   F  #faces by 3 list of mesh face indices
+  // Outputs:
+  //   G  #faces*dim by #V Gradient operator
+  //
+  
+  // Gradient of a scalar function defined on piecewise linear elements (mesh)
+  // is constant on each triangle i,j,k:
+  // grad(Xijk) = (Xj-Xi) * (Vi - Vk)^R90 / 2A + (Xk-Xi) * (Vj - Vi)^R90 / 2A
+  // where Xi is the scalar value at vertex i, Vi is the 3D position of vertex
+  // i, and A is the area of triangle (i,j,k). ^R90 represent a rotation of 
+  // 90 degrees
+  //
+  template <typename T, typename S>
+  IGL_INLINE void gradMat(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &V,
+    const Eigen::Matrix<S, Eigen::Dynamic, Eigen::Dynamic> &F,
+  Eigen::SparseMatrix<T> &G);
+}
+
+#ifdef IGL_HEADER_ONLY
+#  include "gradMat.cpp"
+#endif
+
+#endif