Implement full type mass matrix (#2193)

include/igl/massmatrix.cpp

@@ -11,6 +11,7 @@
 #include "normalize_row_sums.h"
 #include "sparse.h"
 #include "doublearea.h"
+#include "volume.h"
 #include "repmat.h"
 #include <Eigen/Geometry>
 #include <iostream>
@@ -36,9 +37,6 @@ IGL_INLINE void igl::massmatrix(
     eff_type = (simplex_size == 3?MASSMATRIX_TYPE_VORONOI:MASSMATRIX_TYPE_BARYCENTRIC);
   }
 
-  // Not yet supported
-  assert(type!=MASSMATRIX_TYPE_FULL);
-
   if(simplex_size == 3)
   {
     // Triangles
@@ -48,30 +46,48 @@ IGL_INLINE void igl::massmatrix(
     return massmatrix_intrinsic(l,F,type,M);
   }else if(simplex_size == 4)
   {
+    // tetrahedra
+    assert(V.cols() == 3 && "vertices must be defined in 3D for tetrahedra");
+
+    Matrix<Scalar,Dynamic,1> vol;
+    volume(V,F,vol);
+    vol = vol.array().abs();
+
     Matrix<typename DerivedF::Scalar,Dynamic,1> MI;
     Matrix<typename DerivedF::Scalar,Dynamic,1> MJ;
     Matrix<Scalar,Dynamic,1> MV;
-    assert(V.cols() == 3);
-    assert(eff_type == MASSMATRIX_TYPE_BARYCENTRIC);
-    MI.resize(m*4,1); MJ.resize(m*4,1); MV.resize(m*4,1);
-    MI.block(0*m,0,m,1) = F.col(0);
-    MI.block(1*m,0,m,1) = F.col(1);
-    MI.block(2*m,0,m,1) = F.col(2);
-    MI.block(3*m,0,m,1) = F.col(3);
-    MJ = MI;
-    // loop over tets
-    for(int i = 0;i<m;i++)
+
+    switch (eff_type)
     {
-      // http://en.wikipedia.org/wiki/Tetrahedron#Volume
-      Matrix<Scalar,3,1> v0m3,v1m3,v2m3;
-      v0m3.head(V.cols()) = V.row(F(i,0)) - V.row(F(i,3));
-      v1m3.head(V.cols()) = V.row(F(i,1)) - V.row(F(i,3));
-      v2m3.head(V.cols()) = V.row(F(i,2)) - V.row(F(i,3));
-      Scalar v = fabs(v0m3.dot(v1m3.cross(v2m3)))/6.0;
-      MV(i+0*m) = v/4.0;
-      MV(i+1*m) = v/4.0;
-      MV(i+2*m) = v/4.0;
-      MV(i+3*m) = v/4.0;
+      case MASSMATRIX_TYPE_BARYCENTRIC:
+        MI.resize(m*4,1); MJ.resize(m*4,1); MV.resize(m*4,1);
+        MI.block(0*m,0,m,1) = F.col(0);
+        MI.block(1*m,0,m,1) = F.col(1);
+        MI.block(2*m,0,m,1) = F.col(2);
+        MI.block(3*m,0,m,1) = F.col(3);
+        MJ = MI;
+        repmat(vol,4,1,MV);
+        assert(MV.rows()==m*4&&MV.cols()==1);
+        MV.array() /= 4.;
+        break;
+      case MASSMATRIX_TYPE_VORONOI:
+        {
+          assert(false && "Implementation incomplete");
+          break;
+        }
+      case MASSMATRIX_TYPE_FULL:
+        MI.resize(m*16,1); MJ.resize(m*16,1); MV.resize(m*16,1);
+        // indicies and values of the element mass matrix entries in the order
+        // (1,0),(2,0),(3,0),(2,1),(3,1),(0,1),(3,2),(0,2),(1,2),(0,3),(1,3),(2,3),(0,0),(1,1),(2,2),(3,3);
+        MI<<F.col(1),F.col(2),F.col(3),F.col(2),F.col(3),F.col(0),F.col(3),F.col(0),F.col(1),F.col(0),F.col(1),F.col(2),F.col(0),F.col(1),F.col(2),F.col(3);
+        MJ<<F.col(0),F.col(0),F.col(0),F.col(1),F.col(1),F.col(1),F.col(2),F.col(2),F.col(2),F.col(3),F.col(3),F.col(3),F.col(0),F.col(1),F.col(2),F.col(3);
+        repmat(vol,16,1,MV);
+        assert(MV.rows()==m*16&&MV.cols()==1);
+        MV.block(0*m,0,12*m,1) /= 20.;
+        MV.block(12*m,0,4*m,1) /= 10.;
+        break;
+      default:
+        assert(false && "Unknown Mass matrix eff_type");
     }
     sparse(MI,MJ,MV,n,n,M);
   }else

include/igl/massmatrix.h

@@ -21,7 +21,7 @@ namespace igl
     MASSMATRIX_TYPE_VORONOI = 1,
     MASSMATRIX_TYPE_FULL = 2,
     MASSMATRIX_TYPE_DEFAULT = 3,
-    NUM_MASSMATRIX_TYPE = 4
+    NUM_MASSMATRIX_TYPES = 4
   };
 
   // Constructs the mass (area) matrix for a given mesh (V,F).
@@ -36,9 +36,9 @@ namespace igl
   //   V  #V by dim list of mesh vertex positions
   //   F  #F by simplex_size list of mesh elements (triangles or tetrahedra)
   //   type  one of the following ints:
-  //     MASSMATRIX_TYPE_BARYCENTRIC  barycentric
-  //     MASSMATRIX_TYPE_VORONOI voronoi-hybrid {default}
-  //     MASSMATRIX_TYPE_FULL full {not implemented}
+  //     MASSMATRIX_TYPE_BARYCENTRIC  barycentric {default for tetrahedra}
+  //     MASSMATRIX_TYPE_VORONOI voronoi-hybrid {default for triangles, not implemented for tetrahedra}
+  //     MASSMATRIX_TYPE_FULL full
   // Outputs: 
   //   M  #V by #V mass matrix
   //

include/igl/massmatrix_intrinsic.cpp

@@ -111,7 +111,14 @@ IGL_INLINE void igl::massmatrix_intrinsic(
         break;
       }
     case MASSMATRIX_TYPE_FULL:
-      assert(false && "Implementation incomplete");
+      MI.resize(m*9,1); MJ.resize(m*9,1); MV.resize(m*9,1);
+      // indicies and values of the element mass matrix entries in the order
+      // (0,1),(1,0),(1,2),(2,1),(2,0),(0,2),(0,0),(1,1),(2,2);
+      MI<<F.col(0),F.col(1),F.col(1),F.col(2),F.col(2),F.col(0),F.col(0),F.col(1),F.col(2);
+      MJ<<F.col(1),F.col(0),F.col(2),F.col(1),F.col(0),F.col(2),F.col(0),F.col(1),F.col(2);
+      repmat(dblA,9,1,MV);
+      MV.block(0*m,0,6*m,1) /= 24.0;
+      MV.block(6*m,0,3*m,1) /= 12.0;
       break;
     default:
       assert(false && "Unknown Mass matrix eff_type");

include/igl/massmatrix_intrinsic.h

@@ -24,7 +24,7 @@ namespace igl
   //   type  one of the following ints:
   //     MASSMATRIX_TYPE_BARYCENTRIC  barycentric
   //     MASSMATRIX_TYPE_VORONOI voronoi-hybrid {default}
-  //     MASSMATRIX_TYPE_FULL full {not implemented}
+  //     MASSMATRIX_TYPE_FULL full
   // Outputs: 
   //   M  #V by #V mass matrix
   //

include/igl/volume.cpp

@@ -123,4 +123,8 @@ template void igl::volume<Eigen::Matrix<double, -1, 6, 0, -1, 6>, Eigen::Matrix<
 template void igl::volume<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::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 template Eigen::Matrix<double, 1, 3, 1, 1, 3>::Scalar igl::volume_single<Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3> >(Eigen::Matrix<double, 1, 3, 1, 1, 3> const&, Eigen::Matrix<double, 1, 3, 1, 1, 3> const&, Eigen::Matrix<double, 1, 3, 1, 1, 3> const&, Eigen::Matrix<double, 1, 3, 1, 1, 3> const&);
 template void igl::volume<Eigen::Matrix<double,-1,3,0,-1,3>,Eigen::Matrix<int,-1,3,0,-1,3>,Eigen::Matrix<double,-1,1,0,-1,1> >(Eigen::MatrixBase<Eigen::Matrix<double,-1,3,0,-1,3> > const &,Eigen::MatrixBase<Eigen::Matrix<int,-1,3,0,-1,3> > const &,Eigen::PlainObjectBase<Eigen::Matrix<double,-1,1,0,-1,1> > &);
+template void igl::volume<class Eigen::Matrix<double, -1, 3, 0, -1, 3>, class Eigen::Matrix<int, -1, -1, 0, -1, -1>, class Eigen::Matrix<double, -1, 1, 0, -1, 1> >(class Eigen::MatrixBase<class Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, class Eigen::MatrixBase<class Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, class Eigen::PlainObjectBase<class Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template void igl::volume<class Eigen::Matrix<double, -1, -1, 1, -1, -1>,class Eigen::Matrix<int, -1, -1, 0, -1, -1>,class Eigen::Matrix<double, -1, 1, 0, -1, 1> >(class Eigen::MatrixBase<class Eigen::Matrix<double, -1, -1, 1, -1, -1> > const &,class Eigen::MatrixBase<class Eigen::Matrix<int, -1, -1, 0, -1, -1> > const &,class Eigen::PlainObjectBase<class Eigen::Matrix<double, -1, 1, 0, -1, 1> > &);
+template void igl::volume<class Eigen::Matrix<double, -1, 3, 0, -1, 3>, class Eigen::Matrix<int, -1, 4, 0, -1, 4>, class Eigen::Matrix<double, -1, 1, 0, -1, 1> >(class Eigen::MatrixBase<class Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, class Eigen::MatrixBase<class Eigen::Matrix<int, -1, 4, 0, -1, 4> > const&, class Eigen::PlainObjectBase<class Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template void igl::volume<class Eigen::Matrix<double, -1, 3, 1, -1, 3>, class Eigen::Matrix<int, -1, 3, 1, -1, 3>, class Eigen::Matrix<double, -1, 1, 0, -1, 1> >(class Eigen::MatrixBase<class Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, class Eigen::MatrixBase<class Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, class Eigen::PlainObjectBase<class Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/volume.h

@@ -18,7 +18,7 @@ namespace igl
   //
   // Inputs:
   //   V  #V by dim list of vertex positions
-  //   T  #V by 4 list of tet indices
+  //   T  #T by 4 list of tet indices
   // Outputs:
   //   vol  #T list of tetrahedron volumes
   //

tests/include/igl/massmatrix.cpp

@@ -0,0 +1,111 @@
+#include <test_common.h>
+#include <igl/massmatrix.h>
+#include <igl/doublearea.h>
+#include <igl/edges.h>
+#include <igl/vertex_triangle_adjacency.h>
+
+TEST_CASE("massmatrix: full", "[igl]" )
+{
+  const auto test_case = [](const std::string &param)
+  {
+    const double epsilon = 1e-15;
+
+    Eigen::MatrixXd V;
+    Eigen::MatrixXi F;
+    // Load example mesh: GetParam() will be name of mesh file
+    igl::read_triangle_mesh(test_common::data_path(param), V, F);
+    Eigen::SparseMatrix<double> M;
+    igl::massmatrix(V,F,igl::MASSMATRIX_TYPE_FULL,M);
+    REQUIRE(M.rows() == V.rows());
+    REQUIRE(M.cols() == V.rows());
+    Eigen::MatrixXd dblA;
+    igl::doublearea(V,F,dblA);
+    REQUIRE(M.sum() == Approx(dblA.sum()/2.).margin(epsilon));
+  };
+
+  test_common::run_test_cases(test_common::all_meshes(), test_case);
+}
+
+TEST_CASE("massmatrix: barycentric", "[igl]" )
+{
+  const auto test_case = [](const std::string &param)
+  {
+    const double epsilon = 1e-15;
+
+    Eigen::MatrixXd V;
+    Eigen::MatrixXi F;
+    // Load example mesh: GetParam() will be name of mesh file
+    igl::read_triangle_mesh(test_common::data_path(param), V, F);
+    Eigen::SparseMatrix<double> M;
+    igl::massmatrix(V,F,igl::MASSMATRIX_TYPE_BARYCENTRIC,M);
+    REQUIRE(M.rows() == V.rows());
+    REQUIRE(M.cols() == V.rows());
+    Eigen::MatrixXd dblA;
+    igl::doublearea(V,F,dblA);
+    REQUIRE(M.sum() == Approx(dblA.sum()/2.).margin(epsilon));
+  };
+
+  test_common::run_test_cases(test_common::all_meshes(), test_case);
+}
+
+TEST_CASE("massmatrix: cube_barycentric", "[igl]")
+{
+  //The allowed error for this test
+  const double epsilon = 1e-15;
+
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  //This is a cube of dimensions 1.0x1.0x1.0
+  igl::read_triangle_mesh(test_common::data_path("cube.obj"), V, F);
+
+  Eigen::SparseMatrix<double> M;
+
+  //Check the mass matrix of the cube
+  igl::massmatrix(V,F,igl::MASSMATRIX_TYPE_BARYCENTRIC,M);
+
+  REQUIRE(M.rows() == V.rows());
+  REQUIRE(M.cols() == V.rows());
+  //All triangles' areas are 0.5
+  //barycentric area for a vertex is {number of adj triangles} * 0.5 / 3
+  Eigen::VectorXi adj(8);
+  adj << 6,4,4,4,4,4,6,4;
+  for(int f = 0;f<M.rows();f++)
+  {
+    REQUIRE(M.coeff(f,f) == Approx(0.5/3*adj(f)).margin(epsilon));
+  }
+}
+
+TEST_CASE("massmatrix: cube_full", "[igl]")
+{
+  //The allowed error for this test
+  const double epsilon = 1e-15;
+
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  //Check the mass matrix of the cube
+  igl::read_triangle_mesh(test_common::data_path("cube.obj"), V, F);
+
+  Eigen::SparseMatrix<double> M;
+
+  //Check the regular tetrahedron of side sqrt(2)
+  igl::massmatrix(V,F,igl::MASSMATRIX_TYPE_FULL,M);
+
+  REQUIRE(M.rows() == V.rows());
+  REQUIRE(M.cols() == V.rows());
+  Eigen::VectorXi adj(8);
+  adj << 6,4,4,4,4,4,6,4;
+  //All triangles' areas are 0.5
+  //full mass matrix on diagnol is {number of adj triangles} * 0.5 / 6
+  for(int f=0;f<M.rows();++f)
+  {
+    REQUIRE(M.coeff(f,f) == Approx(0.5/6*adj(f)).margin(epsilon));
+  }
+  for (int i=0;i<F.rows();++i)
+  {
+    for (int j=0;j<3;++j)
+    {
+      REQUIRE(M.coeff(F(i,j),F(i,(j+1)%3)) == Approx(0.5/6).margin(epsilon));
+      REQUIRE(M.coeff(F(i,(j+1)%3),F(i,j)) == Approx(0.5/6).margin(epsilon));
+    }
+  }
+}