consolidated igl::grad

Former-commit-id: c8155f94dd2c47a31463ecac60fe107c33d53ca9

include/igl/grad.cpp

@@ -10,14 +10,116 @@
 #include <vector>
 
 #include "per_face_normals.h"
+#include "volume.h"
+#include "doublearea.h"
 
 template <typename DerivedV, typename DerivedF>
-IGL_INLINE void igl::grad(const Eigen::PlainObjectBase<DerivedV>&V,
+IGL_INLINE void grad_tet(const Eigen::PlainObjectBase<DerivedV>&V,
+                     const Eigen::PlainObjectBase<DerivedF>&T,
+                            Eigen::SparseMatrix<typename DerivedV::Scalar> &G,
+                            bool uniform) {
+  using namespace Eigen;
+  assert(T.cols() == 4);
+  const int n = V.rows(); int m = T.rows();
+
+  /*
+      F = [ ...
+      T(:,1) T(:,2) T(:,3); ...
+      T(:,1) T(:,3) T(:,4); ...
+      T(:,1) T(:,4) T(:,2); ...
+      T(:,2) T(:,4) T(:,3)]; */
+  MatrixXi F(4*m,3);
+  for (int i = 0; i < m; i++) {
+    F.row(0*m + i) << T(i,0), T(i,1), T(i,2);
+    F.row(1*m + i) << T(i,0), T(i,2), T(i,3);
+    F.row(2*m + i) << T(i,0), T(i,3), T(i,1);
+    F.row(3*m + i) << T(i,1), T(i,3), T(i,2);
+  }
+  // compute volume of each tet
+  VectorXd vol; igl::volume(V,T,vol);
+
+  VectorXd A(F.rows());
+  MatrixXd N(F.rows(),3);
+  if (!uniform) {
+    // compute tetrahedron face normals
+    igl::per_face_normals(V,F,N); int norm_rows = N.rows();
+    for (int i = 0; i < norm_rows; i++)
+      N.row(i) /= N.row(i).norm();
+    igl::doublearea(V,F,A); A/=2.;
+  } else {
+    // Use a uniform tetrahedra as a reference, with the same volume as the original one:
+    //
+    // Use normals of the uniform tet (V = h*[0,0,0;1,0,0;0.5,sqrt(3)/2.,0;0.5,sqrt(3)/6.,sqrt(2)/sqrt(3)])
+    //         0         0    1.0000
+    //         0.8165   -0.4714   -0.3333
+    //         0          0.9428   -0.3333
+    //         -0.8165   -0.4714   -0.3333
+    for (int i = 0; i < m; i++) {
+      N.row(0*m+i) << 0,0,1;
+      double a = sqrt(2)*std::cbrt(3*vol(i)); // area of a face in a uniform tet with volume = vol(i)
+      A(0*m+i) = (pow(a,2)*sqrt(3))/4.;
+    }
+    for (int i = 0; i < m; i++) {
+      N.row(1*m+i) << 0.8165,-0.4714,-0.3333;
+      double a = sqrt(2)*std::cbrt(3*vol(i));
+      A(1*m+i) = (pow(a,2)*sqrt(3))/4.;
+    }
+    for (int i = 0; i < m; i++) {
+      N.row(2*m+i) << 0,0.9428,-0.3333;
+      double a = sqrt(2)*std::cbrt(3*vol(i));
+      A(2*m+i) = (pow(a,2)*sqrt(3))/4.;
+    }
+    for (int i = 0; i < m; i++) {
+      N.row(3*m+i) << -0.8165,-0.4714,-0.3333;
+      double a = sqrt(2)*std::cbrt(3*vol(i));
+      A(3*m+i) = (pow(a,2)*sqrt(3))/4.;
+    }
+
+  }
+
+  /*  G = sparse( ...
+      [0*m + repmat(1:m,1,4) ...
+       1*m + repmat(1:m,1,4) ...
+       2*m + repmat(1:m,1,4)], ...
+      repmat([T(:,4);T(:,2);T(:,3);T(:,1)],3,1), ...
+      repmat(A./(3*repmat(vol,4,1)),3,1).*N(:), ...
+      3*m,n);*/
+  std::vector<Triplet<double> > G_t;
+  for (int i = 0; i < 4*m; i++) {
+    int T_j; // j indexes : repmat([T(:,4);T(:,2);T(:,3);T(:,1)],3,1)
+    switch (i/m) {
+      case 0:
+        T_j = 3;
+        break;
+      case 1:
+        T_j = 1;
+        break;
+      case 2:
+        T_j = 2;
+        break;
+      case 3:
+        T_j = 0;
+        break;
+    }
+    int i_idx = i%m;
+    int j_idx = T(i_idx,T_j);
+
+    double val_before_n = A(i)/(3*vol(i_idx));
+    G_t.push_back(Triplet<double>(0*m+i_idx, j_idx, val_before_n * N(i,0)));
+    G_t.push_back(Triplet<double>(1*m+i_idx, j_idx, val_before_n * N(i,1)));
+    G_t.push_back(Triplet<double>(2*m+i_idx, j_idx, val_before_n * N(i,2)));
+  }
+  G.resize(3*m,n);
+  G.setFromTriplets(G_t.begin(), G_t.end());
+}
+
+template <typename DerivedV, typename DerivedF>
+IGL_INLINE void grad_tri(const Eigen::PlainObjectBase<DerivedV>&V,
                      const Eigen::PlainObjectBase<DerivedF>&F,
                     Eigen::SparseMatrix<typename DerivedV::Scalar> &G,
                     bool uniform)
 {
-  Eigen::Matrix<typename DerivedV::Scalar,Eigen::Dynamic,3> 
+  Eigen::Matrix<typename DerivedV::Scalar,Eigen::Dynamic,3>
     eperp21(F.rows(),3), eperp13(F.rows(),3);
 
   for (int i=0;i<F.rows();++i)
@@ -51,7 +153,7 @@ IGL_INLINE void igl::grad(const Eigen::PlainObjectBase<DerivedV>&V,
       v1 << 0,0,0;
       v2 << h,0,0;
       v3 << h/2.,(sqrt(3)/2.)*h,0;
-      
+
       // now fix v32,v13,v21 and the normal
       v32 = v3-v2;
       v13 = v1-v3;
@@ -120,109 +222,20 @@ IGL_INLINE void igl::grad(const Eigen::PlainObjectBase<DerivedV>&V,
 }
 
 template <typename DerivedV, typename DerivedF>
-IGL_INLINE void igl::grad_tet(const Eigen::PlainObjectBase<DerivedV>&V,
-                     const Eigen::PlainObjectBase<DerivedF>&T,
-                            Eigen::SparseMatrix<typename DerivedV::Scalar> &G,
-                            bool uniform) {
-  using namespace Eigen;
-  assert(T.cols() == 4);
-  const int n = V.rows(); int m = T.rows();
-
-  /*
-      F = [ ...
-      T(:,1) T(:,2) T(:,3); ...
-      T(:,1) T(:,3) T(:,4); ...
-      T(:,1) T(:,4) T(:,2); ...
-      T(:,2) T(:,4) T(:,3)]; */
-  MatrixXi F(4*m,3);
-  for (int i = 0; i < m; i++) {
-    F.row(0*m + i) << T(i,0), T(i,1), T(i,2);
-    F.row(1*m + i) << T(i,0), T(i,2), T(i,3);
-    F.row(2*m + i) << T(i,0), T(i,3), T(i,1);
-    F.row(3*m + i) << T(i,1), T(i,3), T(i,2);
-  }
-  // compute volume of each tet
-  VectorXd vol; igl::volume(V,T,vol);
-
-  VectorXd A(F.rows());
-  MatrixXd N(F.rows(),3);
-  if (!uniform) {
-    // compute tetrahedron face normals
-    igl::per_face_normals(V,F,N); int norm_rows = N.rows();
-    for (int i = 0; i < norm_rows; i++)
-      N.row(i) /= N.row(i).norm();
-    igl::doublearea(V,F,A); A/=2.;
-  } else {
-    // Use a uniform tetrahedra as a reference, with the same volume as the original one:
-    //
-    // Use normals of the uniform tet (V = h*[0,0,0;1,0,0;0.5,sqrt(3)/2.,0;0.5,sqrt(3)/6.,sqrt(2)/sqrt(3)])
-    //         0         0    1.0000
-    //         0.8165   -0.4714   -0.3333
-    //         0          0.9428   -0.3333
-    //         -0.8165   -0.4714   -0.3333
-    for (int i = 0; i < m; i++) {
-      N.row(0*m+i) << 0,0,1;
-      double a = sqrt(2)*std::cbrt(3*vol(i)); // area of a face in a uniform tet with volume = vol(i)
-      A(0*m+i) = (pow(a,2)*sqrt(3))/4.;
-    }
-    for (int i = 0; i < m; i++) {
-      N.row(1*m+i) << 0.8165,-0.4714,-0.3333;
-      double a = sqrt(2)*std::cbrt(3*vol(i));
-      A(1*m+i) = (pow(a,2)*sqrt(3))/4.;
-    }
-    for (int i = 0; i < m; i++) {
-      N.row(2*m+i) << 0,0.9428,-0.3333;
-      double a = sqrt(2)*std::cbrt(3*vol(i));
-      A(2*m+i) = (pow(a,2)*sqrt(3))/4.;
-    }
-    for (int i = 0; i < m; i++) {
-      N.row(3*m+i) << -0.8165,-0.4714,-0.3333;
-      double a = sqrt(2)*std::cbrt(3*vol(i));
-      A(3*m+i) = (pow(a,2)*sqrt(3))/4.;
-    }
-
-  }
-
-  /*  G = sparse( ...
-      [0*m + repmat(1:m,1,4) ...
-       1*m + repmat(1:m,1,4) ...
-       2*m + repmat(1:m,1,4)], ...
-      repmat([T(:,4);T(:,2);T(:,3);T(:,1)],3,1), ...
-      repmat(A./(3*repmat(vol,4,1)),3,1).*N(:), ...
-      3*m,n);*/
-  std::vector<Triplet<double> > G_t;
-  for (int i = 0; i < 4*m; i++) {
-    int T_j; // j indexes : repmat([T(:,4);T(:,2);T(:,3);T(:,1)],3,1)
-    switch (i/m) {
-      case 0:
-        T_j = 3;
-        break;
-      case 1:
-        T_j = 1;
-        break;
-      case 2:
-        T_j = 2;
-        break;
-      case 3:
-        T_j = 0;
-        break;
-    }
-    int i_idx = i%m;
-    int j_idx = T(i_idx,T_j);
-
-    double val_before_n = A(i)/(3*vol(i_idx));
-    G_t.push_back(Triplet<double>(0*m+i_idx, j_idx, val_before_n * N(i,0)));
-    G_t.push_back(Triplet<double>(1*m+i_idx, j_idx, val_before_n * N(i,1)));
-    G_t.push_back(Triplet<double>(2*m+i_idx, j_idx, val_before_n * N(i,2)));
-  }
-  G.resize(3*m,n);
-  G.setFromTriplets(G_t.begin(), G_t.end());
+IGL_INLINE void igl::grad(const Eigen::PlainObjectBase<DerivedV>&V,
+                     const Eigen::PlainObjectBase<DerivedF>&F,
+                    Eigen::SparseMatrix<typename DerivedV::Scalar> &G,
+                    bool uniform)
+{
+  assert(F.cols() == 3 || F.cols() == 4);
+  if (F.cols() == 3)
+    return grad_tri(V,F,G,uniform);
+  if (F.cols() == 4)
+    return grad_tet(V,F,G,uniform);
 }
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
-// template void igl::grad<double, int>(Eigen::Matrix<double, -1, -1, 0, -1,-1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&,Eigen::SparseMatrix<double, 0, int>&);
-template void igl::grad<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::SparseMatrix<Eigen::Matrix<double, -1, 3, 0, -1, 3>::Scalar, 0, int>&);
-//template void igl::grad<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::SparseMatrix<Eigen::Matrix<double, -1, 3, 0, -1, 3>::Scalar, 0, int>&);
-template void igl::grad<Eigen::Matrix<double, -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&, Eigen::SparseMatrix<Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, 0, int>&);
+template void igl::grad<Eigen::Matrix<double, -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&, Eigen::SparseMatrix<Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, 0, int>&, bool);
+template void igl::grad<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::SparseMatrix<Eigen::Matrix<double, -1, 3, 0, -1, 3>::Scalar, 0, int>&, bool);
 #endif

include/igl/grad.h

@@ -20,14 +20,14 @@ namespace igl {
   //
   // Inputs:
   //   V          #vertices by 3 list of mesh vertex positions
-  //   F          #faces by 3 list of mesh face indices
+  //   F          #faces by 3 list of mesh face indices [or a #faces by 4 list of tetrahedral indices]
   //   uniform    boolean (default false) - Use a uniform mesh instead of the vertices V
   // 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:
+  // is constant on each triangle [tetrahedron] 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
@@ -38,27 +38,7 @@ IGL_INLINE void grad(const Eigen::PlainObjectBase<DerivedV>&V,
                      const Eigen::PlainObjectBase<DerivedF>&F,
                     Eigen::SparseMatrix<typename DerivedV::Scalar> &G,
                     bool uniform = false);
-
-
- // G = grad(V,F)
-  //
-  // Compute the numerical gradient operator for a tet mesh
-  //
-  // Inputs:
-  //   V        #vertices by 3 list of mesh vertex positions
-  //   T        #tets by 4 list of tet indices
-  //   uniform  boolean (default false) - Use a uniform mesh instead of the vertices V
-  // Outputs:
-  //   G  #faces*dim by #V Gradient operator
-  //
-template <typename DerivedV, typename DerivedF>
-IGL_INLINE void grad_tet(const Eigen::PlainObjectBase<DerivedV>&V,
-                     const Eigen::PlainObjectBase<DerivedF>&T,
-                            Eigen::SparseMatrix<typename DerivedV::Scalar> &G,
-                            bool uniform = false);
-
 }
-
 #ifndef IGL_STATIC_LIBRARY
 #  include "grad.cpp"
 #endif