hessian.cpp 2.1 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455
  1. // This file is part of libigl, a simple c++ geometry processing library.
  2. //
  3. // Copyright (C) 2017 Alec Jacobson <alecjacobson@gmail.com> and Oded Stein <oded.stein@columbia.edu>
  4. //
  5. // This Source Code Form is subject to the terms of the Mozilla Public License
  6. // v. 2.0. If a copy of the MPL was not distributed with this file, You can
  7. // obtain one at http://mozilla.org/MPL/2.0/.
  8. #include "hessian.h"
  9. #include <vector>
  10. #include <igl/grad.h>
  11. #include <igl/doublearea.h>
  12. #include <igl/repdiag.h>
  13. template <typename DerivedV, typename DerivedF, typename Scalar>
  14. IGL_INLINE void igl::hessian(
  15. const Eigen::MatrixBase<DerivedV> & V,
  16. const Eigen::MatrixBase<DerivedF> & F,
  17. Eigen::SparseMatrix<Scalar>& H)
  18. {
  19. typedef typename DerivedV::Scalar denseScalar;
  20. typedef typename Eigen::Matrix<denseScalar, Eigen::Dynamic, 1> VecXd;
  21. typedef typename Eigen::SparseMatrix<Scalar> SparseMat;
  22. typedef typename Eigen::DiagonalMatrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> DiagMat;
  23. int dim = V.cols();
  24. assert((dim==2 || dim==3) && "The dimension of the vertices should be 2 or 3");
  25. //Construct the combined gradient matric
  26. SparseMat G;
  27. igl::grad(Eigen::PlainObjectBase<DerivedV>(V), Eigen::PlainObjectBase<DerivedF>(F), G, false);
  28. SparseMat GG(F.rows(), dim*V.rows());
  29. GG.reserve(G.nonZeros());
  30. for(int i=0; i<dim; ++i)
  31. GG.middleCols(i*G.cols(), G.cols()) = G.middleRows(i*F.rows(), F.rows());
  32. SparseMat D;
  33. igl::repdiag(GG,dim,D);
  34. //Compute area matrix
  35. VecXd areas;
  36. igl::doublearea(V, F, areas);
  37. DiagMat A = (0.5*areas).replicate(dim,1).asDiagonal();
  38. //Compute FEM Hessian
  39. H = D.transpose()*A*G;
  40. }
  41. #ifdef IGL_STATIC_LIBRARY
  42. // Explicit template instantiation
  43. template void igl::hessian<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&, Eigen::SparseMatrix<double, 0, int>&);
  44. #endif