hybrid mass matrix for tets (#2364)

Co-authored-by: Alec Jacobson <[email protected]>

include/igl/centroid.cpp

@@ -70,6 +70,8 @@ IGL_INLINE void igl::centroid(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::centroid<Eigen::Matrix<double, 8, 3, 0, 8, 3>, Eigen::Matrix<int, 12, 3, 0, 12, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, double>(Eigen::MatrixBase<Eigen::Matrix<double, 8, 3, 0, 8, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, 12, 3, 0, 12, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&, double&);
 template void igl::centroid<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, double>(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, 3, 1, 1, 3> >&, double&);
 // generated by autoexplicit.sh
 template void igl::centroid<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3>, Eigen::Matrix<float, 3, 1, 0, 3, 1>, float>(Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, 3, 1, 0, 3, 1> >&, float&);

include/igl/circumradius.cpp

@@ -8,6 +8,7 @@
 #include "circumradius.h"
 #include "edge_lengths.h"
 #include "doublearea.h"
+#include <Eigen/QR>
 template <
   typename DerivedV, 
   typename DerivedF,
@@ -17,6 +18,9 @@ IGL_INLINE void igl::circumradius(
   const Eigen::MatrixBase<DerivedF> & F,
   Eigen::PlainObjectBase<DerivedR> & R)
 {
+  // cols at compiletime should be Dynamic or 3
+  static_assert(DerivedF::ColsAtCompileTime == Eigen::Dynamic || DerivedF::ColsAtCompileTime == 3,"F should contain triangles");
+  assert(F.cols() == 3 && "F should contain triangles");
   Eigen::Matrix<typename DerivedV::Scalar,Eigen::Dynamic,3> l;
   igl::edge_lengths(V,F,l);
   DerivedR A;
@@ -25,6 +29,59 @@ IGL_INLINE void igl::circumradius(
   R = l.col(0).array() * l.col(1).array() * l.col(2).array() / (2.0*A.array());
 }
 
+template <
+  typename DerivedV, 
+  typename DerivedT,
+  typename DerivedR,
+  typename DerivedC,
+  typename DerivedB>
+IGL_INLINE void igl::circumradius(
+  const Eigen::MatrixBase<DerivedV> & V, 
+  const Eigen::MatrixBase<DerivedT> & T,
+  Eigen::PlainObjectBase<DerivedR> & R,
+  Eigen::PlainObjectBase<DerivedC> & C,
+  Eigen::PlainObjectBase<DerivedB> & B)
+{
+  using Scalar = typename DerivedV::Scalar;
+  const int TCOLS = DerivedT::ColsAtCompileTime;
+  const int VCOLS = DerivedV::ColsAtCompileTime;
+  const int ACOLS = TCOLS==Eigen::Dynamic?Eigen::Dynamic:TCOLS+1;
+  const int ss = T.cols();
+  R.resize(T.rows(),1);
+  C.resize(T.rows(),V.cols());
+  B.resize(T.rows(),T.cols());
+  for(int i = 0;i<T.rows();i++)
+  {
+    Eigen::Matrix<Scalar,ACOLS,ACOLS> A(T.cols()+1,T.cols()+1);
+    // Not sure if this .eval() is a good idea
+    const auto Vi = V(T.row(i),Eigen::all).eval();
+    A.topLeftCorner(T.cols(),T.cols()) = 2*(Vi*Vi.transpose());
+    A.block(0,T.cols(),T.cols(),1).setConstant(1);
+    A.block(T.cols(),0,1,T.cols()).setConstant(1);
+    A(T.cols(),T.cols()) = 0;
+    Eigen::Matrix<Scalar,ACOLS,1> b(T.cols()+1,1);
+    b.head(T.cols()) = (Vi.array().square().rowwise().sum()).eval();
+    b(T.cols()) = 1;
+    // [B;λ] = A\b
+    const auto x = (A.colPivHouseholderQr().solve(b)).eval();
+    B.row(i) = x.head(T.cols());
+    C.row(i) = B.row(i) * Vi;
+    const Scalar lambda = x(T.cols());
+    R(i) = std::sqrt(lambda + C.row(i).squaredNorm());
+  }
+}
+
 #ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::circumradius<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<double, -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> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::circumradius<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 void igl::circumradius<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 4, 0, -1, 4> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 4, 0, -1, 4> >&);
+template void igl::circumradius<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
+template void igl::circumradius<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<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 4, 0, -1, 4> >(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> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 4, 0, -1, 4> >&);
+template void igl::circumradius<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
+template void igl::circumradius<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 4, 0, -1, 4> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 4, 0, -1, 4> >&);
+template void igl::circumradius<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 4, 0, -1, 4>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 4, 0, -1, 4> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 4, 0, -1, 4> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 4, 0, -1, 4> >&);
+template void igl::circumradius<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::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(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> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
 template void igl::circumradius<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> >&);
 #endif

include/igl/circumradius.h

@@ -25,6 +25,26 @@ namespace igl
     const Eigen::MatrixBase<DerivedV> & V, 
     const Eigen::MatrixBase<DerivedF> & F,
     Eigen::PlainObjectBase<DerivedR> & R);
+  /// Generic version 
+  ///
+  /// @param[in] V  #V by dim list of mesh vertex positions
+  /// @param[in] T  #T by simplex-size list of simplex indices into V
+  /// @param[out] R  #T list of circumradius
+  /// @param[out] C  #T by dim list of circumcenter
+  /// @param[out] B  #T by simplex-size list of barycentric coordinates of circumcenter
+  template <
+    typename DerivedV, 
+    typename DerivedT,
+    typename DerivedR,
+    typename DerivedC,
+    typename DerivedB>
+  IGL_INLINE void circumradius(
+    const Eigen::MatrixBase<DerivedV> & V, 
+    const Eigen::MatrixBase<DerivedT> & T,
+    Eigen::PlainObjectBase<DerivedR> & R,
+    Eigen::PlainObjectBase<DerivedC> & C,
+    Eigen::PlainObjectBase<DerivedB> & B);
+
 }
 #ifndef IGL_STATIC_LIBRARY
 #  include "circumradius.cpp"

include/igl/massmatrix.cpp

@@ -11,6 +11,7 @@
 #include "sparse.h"
 #include "doublearea.h"
 #include "volume.h"
+#include "voronoi_mass.h"
 #include "repmat.h"
 #include <Eigen/Geometry>
 #include <iostream>
@@ -71,7 +72,9 @@ IGL_INLINE void igl::massmatrix(
         break;
       case MASSMATRIX_TYPE_VORONOI:
         {
-          assert(false && "Implementation incomplete");
+          MI = decltype(MI)::LinSpaced(n,0,n-1);
+          MJ = MI;
+          voronoi_mass(V,F,MV);
           break;
         }
       case MASSMATRIX_TYPE_FULL:

include/igl/massmatrix.h

@@ -42,7 +42,7 @@ namespace igl
   /// @param[in] F  #F by simplex_size list of mesh elements (triangles or tetrahedra)
   /// @param[in] type  one of the following ints:
   ///     MASSMATRIX_TYPE_BARYCENTRIC  barycentric {default for tetrahedra}
-  ///     MASSMATRIX_TYPE_VORONOI voronoi-hybrid {default for triangles, not implemented for tetrahedra}
+  ///     MASSMATRIX_TYPE_VORONOI voronoi-hybrid {default for triangles}
   ///     MASSMATRIX_TYPE_FULL full
   /// @param[out] M  #V by #V mass matrix
   ///

include/igl/voronoi_mass.cpp

@@ -0,0 +1,142 @@
+#include "voronoi_mass.h"
+#include "circumradius.h"
+#include "centroid.h"
+#include "unique_simplices.h"
+
+template <
+  typename DerivedV,
+  typename DerivedT,
+  typename DerivedM>
+void igl::voronoi_mass(
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedT> & T,
+    Eigen::PlainObjectBase<DerivedM> & M)
+{
+  // DerivedM::Scalar must be same as DerivedV::Scalar
+  static_assert(std::is_same<typename DerivedM::Scalar, typename DerivedV::Scalar>::value, "DerivedM::Scalar must be same as DerivedV::Scalar");
+  // DerivedM must be a vector (rows or cols at compile time = 1)
+  static_assert(DerivedM::RowsAtCompileTime == 1 || DerivedM::ColsAtCompileTime == 1, "DerivedM must be a vector (rows or cols at compile time = 1)");
+  // DerivedT must have Dynamic or 4 cols
+  static_assert(DerivedT::ColsAtCompileTime == Eigen::Dynamic || DerivedT::ColsAtCompileTime == 4, "DerivedT must have Dynamic or 4 cols");
+  assert(T.cols() == 4 && "Tetrahedra should have 4 vertices");
+  // DerivedV must have Dynamic or 3 cols
+  static_assert(DerivedV::ColsAtCompileTime == Eigen::Dynamic || DerivedV::ColsAtCompileTime == 3, "DerivedV must have Dynamic or 3 cols");
+  assert(V.cols() == 3 && "V should have 3 columns");
+
+
+  using Scalar = typename DerivedV::Scalar;
+  using VectorXS = Eigen::Matrix<Scalar,Eigen::Dynamic,1>;
+  using MatrixX3S = Eigen::Matrix<Scalar,Eigen::Dynamic,3>;
+  using MatrixX4S = Eigen::Matrix<Scalar,Eigen::Dynamic,4>;
+  using MatrixX3I = Eigen::Matrix<int,Eigen::Dynamic,3>;
+  using MatrixX4I = Eigen::Matrix<int,Eigen::Dynamic,4>;
+  MatrixX3I F;
+  MatrixX4I I;
+  {
+    MatrixX3I allF(T.rows()*T.cols(),3);
+    for(int i = 0;i<T.rows();i++)
+    {
+      for(int j = 0;j<T.cols();j++)
+      {
+        allF.row(j*T.rows()+i) <<
+          T(i,(j+1)%T.cols()),
+          T(i,(j+2)%T.cols()),
+          T(i,(j+3)%T.cols());
+      }
+    }
+    Eigen::VectorXi _;
+    Eigen::VectorXi Ivec;
+    igl::unique_simplices(allF,F,_,Ivec);
+    I = Ivec.reshaped(T.rows(),T.cols());
+  }
+
+  // Face circumcenters
+  MatrixX3S CF;
+  {
+    VectorXS R;
+    MatrixX3S B;
+    circumradius(V,F,R,CF,B);
+    for(int i = 0;i<F.rows();i++)
+    {
+      int j;
+      // Snap to edge "circumcenter"
+      if(B.row(i).minCoeff(&j) < 0)
+      {
+        CF.row(i) = (V.row(F(i,(j+1)%3)) + V.row(F(i,(j+2)%3)))*0.5;
+      }
+    }
+  }
+  MatrixX3S CT;
+  {
+    VectorXS R;
+    MatrixX4S B;
+    circumradius(V,T,R,CT,B);
+    for(int i = 0;i<T.rows();i++)
+    {
+      int j;
+      // Snap to face "circumcenter"
+      if(B.row(i).minCoeff(&j) < 0)
+      {
+        CT.row(i) = CF.row(I(i,j));
+      }
+    }
+  }
+
+  M.setZero(V.rows());
+  for(int i = 0;i<T.rows();i++)
+  {
+    for(int j = 0;j<4;j++)
+    {
+      Eigen::Matrix<Scalar,8,3> U(8,V.cols());
+      U.row(0) = V.row(T(i,j));
+      // edge circumcenters
+      for(int k = 1;k<4;k++)
+      {
+        U.row(k) = 0.5*(U.row(0) + V.row(T(i,(j+k)%4)));
+      }
+      // face circumcenters
+      {
+        U.block(4,0,3,V.cols()) <<
+          CF.row(I(i,(j+1)%4)),
+          CF.row(I(i,(j+2)%4)),
+          CF.row(I(i,(j+3)%4));
+      }
+      // Tet circumcenter
+      {
+        U.row(7) = CT.row(i);
+      }
+      Eigen::Matrix<int,12,3> Fij(12,3);
+      Fij<< 
+        4,2,0,
+        5,3,0,
+        6,1,0,
+        7,2,4,
+        7,3,5,
+        7,1,6,
+        4,0,3,
+        5,0,1,
+        6,0,2,
+        7,4,3,
+        7,5,1,
+        7,6,2;
+      if(j%2)
+      {
+        Fij = Fij.rowwise().reverse().eval();
+      }
+      Scalar vol;
+      {
+        Eigen::Matrix<Scalar,1,3> cen;
+        igl::centroid(U,Fij,cen,vol);
+      }
+      M(T(i,j)) += vol;
+    }
+  }
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template void igl::voronoi_mass<Eigen::Matrix<double, -1, -1, 1, -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, 1, -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 void igl::voronoi_mass<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 void igl::voronoi_mass<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, 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, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template void igl::voronoi_mass<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 4, 0, -1, 4>, 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, 4, 0, -1, 4> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+#endif

include/igl/voronoi_mass.h

@@ -0,0 +1,37 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2024 Alec Jacobson <[email protected]>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_VORONOI_MASS_H
+#define IGL_VORONOI_MASS_H
+
+#include <igl/igl_inline.h>
+#include <Eigen/Core>
+namespace igl
+{
+/// Compute the mass matrix entries for a given tetrahedral mesh (V,T) using the
+/// "hybrid" voronoi volume of each vertex.
+/// 
+///  @param[in] V  #V by 3 list of vertex positions
+///  @param[in] T  #T by 4 list of element indices into V
+///  @param[out] M  #V list of mass matrix diagonal entries (will be positive as
+///    long as tets are not degenerate)
+///
+template <
+  typename DerivedV,
+  typename DerivedT,
+  typename DerivedM>
+IGL_INLINE void voronoi_mass(
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedT> & T,
+    Eigen::PlainObjectBase<DerivedM> & M);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "voronoi_mass.cpp"
+#endif
+
+#endif

tests/include/igl/circumradius.cpp

@@ -0,0 +1,165 @@
+#include <test_common.h>
+#include <igl/circumradius.h>
+
+TEST_CASE("circumradius: equilateral-triangle", "[igl]" )
+{
+  // Equilateral triangle
+  Eigen::MatrixXd V(3,2);
+  V << 0,0, 1,0, 0.5,sqrt(3.0)/2.0;
+  Eigen::MatrixXi F(1,3);
+  F << 0,1,2;
+  Eigen::VectorXd R;
+  igl::circumradius(V,F,R);
+  REQUIRE (R.size() == 1);
+  REQUIRE (R(0) == Approx(sqrt(3.0)/3.0).margin(1e-15));
+  Eigen::MatrixXd C;
+  Eigen::MatrixXd B;
+  igl::circumradius(V,F,R,C,B);
+  REQUIRE (R.size() == 1);
+  REQUIRE (R(0) == Approx(sqrt(3.0)/3.0).margin(1e-15));
+  REQUIRE (C.rows() == 1);
+  REQUIRE (C.cols() == 2);
+  REQUIRE (C(0,0) == Approx(0.5).margin(1e-15));
+  REQUIRE (C(0,1) == Approx(sqrt(3.0)/6.0).margin(1e-15));
+  REQUIRE (B.rows() == 1);
+  REQUIRE (B.cols() == 3);
+  REQUIRE (B(0,0) == Approx(1.0/3.0).margin(1e-15));
+  REQUIRE (B(0,1) == Approx(1.0/3.0).margin(1e-15));
+  REQUIRE (B(0,2) == Approx(1.0/3.0).margin(1e-15));
+}
+
+TEST_CASE("circumradius: right-triangle", "[igl]" )
+{
+  // Right triangle
+  Eigen::MatrixXd V(3,2);
+  V << 0,0, 1,0, 0,1;
+  Eigen::MatrixXi F(1,3);
+  F << 0,1,2;
+  Eigen::VectorXd R;
+  igl::circumradius(V,F,R);
+  REQUIRE (R.size() == 1);
+  REQUIRE (R(0) == Approx(sqrt(2.0)/2.0).margin(1e-15));
+  Eigen::MatrixXd C;
+  Eigen::MatrixXd B;
+  igl::circumradius(V,F,R,C,B);
+  REQUIRE (R.size() == 1);
+  REQUIRE (R(0) == Approx(sqrt(2.0)/2.0).margin(1e-15));
+  REQUIRE (C.rows() == 1);
+  REQUIRE (C.cols() == 2);
+  REQUIRE (C(0,0) == Approx(0.5).margin(1e-15));
+  REQUIRE (C(0,1) == Approx(0.5).margin(1e-15));
+  REQUIRE (B.rows() == 1);
+  REQUIRE (B.cols() == 3);
+  REQUIRE (B(0,0) == Approx(0.0).margin(1e-15));
+  REQUIRE (B(0,1) == Approx(0.5).margin(1e-15));
+  REQUIRE (B(0,2) == Approx(0.5).margin(1e-15));
+}
+
+// Obtuse
+TEST_CASE("circumradius: obtuse-triangle", "[igl]" )
+{
+  Eigen::MatrixXd V(3,2);
+  V << 0,0, 4,0, 2,1;
+  Eigen::MatrixXi F(1,3);
+  F << 0,1,2;
+  Eigen::VectorXd R;
+  igl::circumradius(V,F,R);
+  REQUIRE (R.size() == 1);
+  REQUIRE (R(0) == Approx(2.5).margin(1e-15));
+
+  Eigen::MatrixXd C;
+  Eigen::MatrixXd B;
+  igl::circumradius(V,F,R,C,B);
+  REQUIRE (R.size() == 1);
+  REQUIRE (R(0) == Approx(2.5).margin(1e-15));
+  REQUIRE (C.rows() == 1);
+  REQUIRE (C.cols() == 2);
+  REQUIRE (C(0,0) == Approx(2).margin(1e-15));
+  REQUIRE (C(0,1) == Approx(-1.5).margin(1e-15));
+  REQUIRE (B.rows() == 1);
+  REQUIRE (B.cols() == 3);
+  REQUIRE (B(0,0) == Approx(1.25).margin(1e-15));
+  REQUIRE (B(0,1) == Approx(1.25).margin(1e-15));
+  REQUIRE (B(0,2) == Approx(-1.5).margin(1e-15));
+}
+
+// Tetrahedra test cases
+
+TEST_CASE("circumradius: equilateral-tetrahedra", "[igl]" )
+{
+  // Tetrahedra
+  Eigen::MatrixXd V(4,3);
+  V << 0,0,0, 1,0,0, 0.5,sqrt(3.0)/2.0,0, 0.5,sqrt(3.0)/6.0,sqrt(2.0/3.0);
+  Eigen::MatrixXi T(1,4);
+  T << 0,1,2,3;
+  Eigen::VectorXd R;
+  Eigen::MatrixXd C;
+  Eigen::MatrixXd B;
+  igl::circumradius(V,T,R,C,B);
+  REQUIRE (R.size() == 1);
+  REQUIRE (R(0) == Approx(sqrt(3.0/8.0)).margin(1e-15));
+  REQUIRE (C.rows() == 1);
+  REQUIRE (C.cols() == 3);
+  REQUIRE (C(0,0) == Approx(0.5).margin(1e-15));
+  REQUIRE (C(0,1) == Approx(sqrt(3.0)/6.0).margin(1e-15));
+  REQUIRE (C(0,2) == Approx(sqrt(1.0/24.0)).margin(1e-15));
+  REQUIRE (B.rows() == 1);
+  REQUIRE (B.cols() == 4);
+  REQUIRE (B(0,0) == Approx(0.25).margin(1e-15));
+  REQUIRE (B(0,1) == Approx(0.25).margin(1e-15));
+  REQUIRE (B(0,2) == Approx(0.25).margin(1e-15));
+  REQUIRE (B(0,3) == Approx(0.25).margin(1e-15));
+}
+
+TEST_CASE("circumradius: right-tetrahedra", "[igl]" )
+{
+  // Tetrahedra
+  Eigen::MatrixXd V(4,3);
+  V << 0,0,0, 1,0,0, 0,1,0, 0,0,1;
+  Eigen::MatrixXi T(1,4);
+  T << 0,1,2,3;
+  Eigen::VectorXd R;
+  Eigen::MatrixXd C;
+  Eigen::MatrixXd B;
+  igl::circumradius(V,T,R,C,B);
+  REQUIRE (R.size() == 1);
+  REQUIRE (R(0) == Approx(sqrt(3.0)/2.0).margin(1e-15));
+  REQUIRE (C.rows() == 1);
+  REQUIRE (C.cols() == 3);
+  REQUIRE (C(0,0) == Approx(0.5).margin(1e-15));
+  REQUIRE (C(0,1) == Approx(0.5).margin(1e-15));
+  REQUIRE (C(0,2) == Approx(0.5).margin(1e-15));
+  REQUIRE (B.rows() == 1);
+  REQUIRE (B.cols() == 4);
+  REQUIRE (B(0,0) == Approx(-0.5).margin(1e-15));
+  REQUIRE (B(0,1) == Approx( 0.5).margin(1e-15));
+  REQUIRE (B(0,2) == Approx( 0.5).margin(1e-15));
+  REQUIRE (B(0,3) == Approx( 0.5).margin(1e-15));
+}
+
+// Obtuse tetrahedron
+
+TEST_CASE("circumradius: obtuse-tetrahedra", "[igl]" )
+{
+  Eigen::MatrixXd V(4,3);
+  V << 0,0,0, 4,0,0, 2,1,0, 2,1,3;
+  Eigen::MatrixXi T(1,4);
+  T << 0,1,2,3;
+  Eigen::VectorXd R;
+  Eigen::MatrixXd C;
+  Eigen::MatrixXd B;
+  igl::circumradius(V,T,R,C,B);
+  REQUIRE (R.size() == 1);
+  REQUIRE (R(0) == Approx(sqrt(17.0/2.0)).margin(1e-15));
+  REQUIRE (C.rows() == 1);
+  REQUIRE (C.cols() == 3);
+  REQUIRE (C(0,0) == Approx(2).margin(1e-15));
+  REQUIRE (C(0,1) == Approx(-1.5).margin(1e-15));
+  REQUIRE (C(0,2) == Approx(1.5).margin(1e-15));
+  REQUIRE (B.rows() == 1);
+  REQUIRE (B.cols() == 4);
+  REQUIRE (B(0,0) == Approx(1.25).margin(1e-15));
+  REQUIRE (B(0,1) == Approx(1.25).margin(1e-15));
+  REQUIRE (B(0,2) == Approx(-2.0).margin(1e-15));
+  REQUIRE (B(0,3) == Approx( 0.5).margin(1e-15));
+}

tests/include/igl/voronoi_mass.cpp

@@ -0,0 +1,60 @@
+#include <test_common.h>
+#include <igl/voronoi_mass.h>
+#include <igl/matlab_format.h>
+#include <iostream>
+
+TEST_CASE("voronoi_mass: equilateral-tetrahedra", "[igl]" )
+{
+  // Tetrahedra
+  Eigen::MatrixXd V(4,3);
+  V << 0,0,0, 1,0,0, 0.5,sqrt(3.0)/2.0,0, 0.5,sqrt(3.0)/6.0,sqrt(2.0/3.0);
+  Eigen::MatrixXi T(1,4);
+  T << 0,1,2,3;
+  Eigen::VectorXd M;
+  igl::voronoi_mass(V,T,M);
+  Eigen::VectorXd Mgt(4);
+  Mgt <<
+  0.0294627825494395,
+  0.0294627825494395,
+  0.0294627825494395,
+  0.0294627825494395;
+  test_common::assert_near(M,Mgt,1e-15);
+}
+
+TEST_CASE("voronoi_mass: right-tetrahedra", "[igl]" )
+{
+  // Tetrahedra
+  Eigen::MatrixXd V(4,3);
+  V << 0,0,0, 1,0,0, 0,1,0, 0,0,1;
+  Eigen::MatrixXi T(1,4);
+  T << 0,1,2,3;
+  Eigen::VectorXd M;
+  igl::voronoi_mass(V,T,M);
+  Eigen::VectorXd Mgt(4);
+  Mgt <<
+  0.0833333333333333,
+  0.0277777777777778,
+  0.0277777777777778,
+  0.0277777777777778;
+  test_common::assert_near(M,Mgt,1e-15);
+}
+
+// Obtuse tetrahedron
+
+TEST_CASE("voronoi_mass: obtuse-tetrahedra", "[igl]" )
+{
+  Eigen::MatrixXd V(4,3);
+  V << 0,0,0, 4,0,0, 2,1,0, 2,1,3;
+  Eigen::MatrixXi T(1,4);
+  T << 0,1,2,3;
+  Eigen::VectorXd M;
+  igl::voronoi_mass(V,T,M);
+  Eigen::VectorXd Mgt(4);
+  Mgt <<
+  0.325,
+  0.325,
+      1,
+   0.35;
+  test_common::assert_near(M,Mgt,1e-15);
+}
+