libigl/include/igl/massmatrix.cpp

// This file is part of libigl, a simple c++ geometry processing library.
// 
// Copyright (C) 2013 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/.
#include "massmatrix.h"
#include "massmatrix_intrinsic.h"
#include "edge_lengths.h"
#include "sparse.h"
#include "doublearea.h"
#include "volume.h"
#include "voronoi_mass.h"
#include "repmat.h"
#include <Eigen/Geometry>
#include <iostream>

// For std::enable_if
#include <type_traits>

namespace igl
{
  template <typename DerivedV, typename DerivedF, typename Scalar, int simplex_size = DerivedF::ColsAtCompileTime>
  struct MassMatrixHelper;
  
  // This would be easier with C++17 if constexpr
  // Specialization for triangles
  template <typename DerivedV, typename DerivedF, typename Scalar>
  struct MassMatrixHelper<DerivedV, DerivedF, Scalar, 3> {
    static void compute(
      const Eigen::MatrixBase<DerivedV>& V, 
      const Eigen::MatrixBase<DerivedF>& F, 
      const MassMatrixType type,
      Eigen::SparseMatrix<Scalar>& M)
    {
      MassMatrixType eff_type = 
        type == MASSMATRIX_TYPE_DEFAULT? MASSMATRIX_TYPE_VORONOI : type;
      Eigen::Matrix<Scalar, Eigen::Dynamic, 3> l;
      igl::edge_lengths(V, F, l);
      return massmatrix_intrinsic(l, F, eff_type, M);
    }
  };
  
  // Specialization for tetrahedra
  template <typename DerivedV, typename DerivedF, typename Scalar>
  struct MassMatrixHelper<DerivedV, DerivedF, Scalar, 4> {
    static void compute(
      const Eigen::MatrixBase<DerivedV>& V, 
      const Eigen::MatrixBase<DerivedF>& F, 
      const MassMatrixType type,
      Eigen::SparseMatrix<Scalar>& M)
    {
      const int n = V.rows();
      const int m = F.rows();
      using Eigen::Matrix;
      using Eigen::Dynamic;
      MassMatrixType eff_type = 
        type == MASSMATRIX_TYPE_DEFAULT? MASSMATRIX_TYPE_BARYCENTRIC: type;
      Eigen::Matrix<Scalar, Eigen::Dynamic, 1> vol;
      volume(V, F, vol);
      vol = vol.array().abs();
      Eigen::Matrix<typename DerivedF::Scalar ,Eigen::Dynamic,1> MI;
      Eigen::Matrix<typename DerivedF::Scalar ,Eigen::Dynamic,1> MJ;
      Eigen::Matrix<Scalar ,Eigen::Dynamic,1> MV;
  
      switch (eff_type)
      {
        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:
          {
            MI = decltype(MI)::LinSpaced(n,0,n-1);
            MJ = MI;
            voronoi_mass(V,F,MV);
            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);
    }
  };
  
  // General template for handling Eigen::Dynamic at runtime
  template <typename DerivedV, typename DerivedF, typename Scalar>
  struct MassMatrixHelper<DerivedV, DerivedF, Scalar, Eigen::Dynamic> {
    static void compute(
      const Eigen::MatrixBase<DerivedV>& V,
      const Eigen::MatrixBase<DerivedF>& F,
      const MassMatrixType type,
      Eigen::SparseMatrix<Scalar>& M)
    {
      if (F.cols() == 3) {
        MassMatrixHelper<DerivedV, DerivedF, Scalar, 3>::compute(V, F, type, M);
      } else if (F.cols() == 4) {
        MassMatrixHelper<DerivedV, DerivedF, Scalar, 4>::compute(V, F, type, M);
      } else {
        // Handle unsupported simplex size at runtime
        assert(false && "Unsupported simplex size");
      }
    }
  };
}

template <typename DerivedV, typename DerivedF, typename Scalar>
IGL_INLINE void igl::massmatrix(
  const Eigen::MatrixBase<DerivedV> & V, 
  const Eigen::MatrixBase<DerivedF> & F, 
  const MassMatrixType type,
  Eigen::SparseMatrix<Scalar>& M)
{
  MassMatrixHelper<DerivedV, DerivedF, Scalar, DerivedF::ColsAtCompileTime>::compute(V, F, type, M);
}

#ifdef IGL_STATIC_LIBRARY
// Explicit template instantiation
// generated by autoexplicit.sh
template void igl::massmatrix<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::MassMatrixType, Eigen::SparseMatrix<double, 0, int>&);
// generated by autoexplicit.sh
template void igl::massmatrix<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::MassMatrixType, Eigen::SparseMatrix<double, 0, int>&);
// generated by autoexplicit.sh
template void igl::massmatrix<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 4, 0, -1, 4>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 4, 0, -1, 4> > const&, igl::MassMatrixType, Eigen::SparseMatrix<double, 0, int>&);
// generated by autoexplicit.sh
template void igl::massmatrix<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, igl::MassMatrixType, Eigen::SparseMatrix<double, 0, int>&);
template void igl::massmatrix<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, igl::MassMatrixType, Eigen::SparseMatrix<double, 0, int>&);
template void igl::massmatrix<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::MassMatrixType, Eigen::SparseMatrix<double, 0, int>&);
#endif