#ifndef IGL_COTMATRIX_H #define IGL_COTMATRIX_H #define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET #include #include // History: // Used const references rather than copying the entire mesh // Alec 9 October 2011 // removed cotan (uniform weights) optional parameter it was building a buggy // half of the uniform laplacian, please see adjacency_matrix istead // Alec 9 October 2011 namespace igl { // Constructs the cotangent stiffness matrix (discrete laplacian) for a given // mesh (V,F). // // Templates: // DerivedV derived type of eigen matrix for V (e.g. derived from // MatirxXd) // DerivedF derived type of eigen matrix for F (e.g. derived from // MatrixXi) // Scalar scalar type for eigen sparse matrix (e.g. double) // Inputs: // V #V by dim list of mesh vertex positions // F #F by simplex_size list of mesh faces (must be triangles) // Outputs: // L #V by #V cotangent matrix, each row i corresponding to V(i,:) // // See also: adjacency_matrix // // Known bugs: off by 1e-16 on regular grid. I think its a problem of // arithmetic order in cotangent.h: C(i,e) = (arithmetic)/dblA/4 template inline void cotmatrix( const Eigen::MatrixBase & V, const Eigen::MatrixBase & F, Eigen::SparseMatrix& L); } // Implementation // For error printing #include #include "cotangent.h" template inline void igl::cotmatrix( const Eigen::MatrixBase & V, const Eigen::MatrixBase & F, Eigen::SparseMatrix& L) { using namespace igl; using namespace Eigen; DynamicSparseMatrix dyn_L (V.rows(), V.rows()); Matrix edges; int simplex_size = F.cols(); // 3 for triangles, 4 for tets assert(simplex_size == 3 || simplex_size == 4); if(simplex_size == 3) { // This is important! it could decrease the comptuation time by a factor of 2 // Laplacian for a closed 2d manifold mesh will have on average 7 entries per // row dyn_L.reserve(7*V.rows()); edges.resize(3,2); edges << 1,2, 2,0, 0,1; }else if(simplex_size == 4) { dyn_L.reserve(17*V.rows()); edges.resize(6,2); edges << 1,2, 2,0, 0,1, 3,0, 3,1, 3,2; }else { return; } // Gather cotangents Matrix C; cotangent(V,F,C); // Loop over triangles for(int i = 0; i < F.rows(); i++) { // loop over edges of element for(int e = 0;e(dyn_L); } #endif