1 // This file is part of libigl, a simple c++ geometry processing library.
2 //
3 // Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com>
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 "LinSpaced.h"
9 #include "normal_derivative.h"
10 #include "cotmatrix_entries.h"
11 #include "slice.h"
12 #include <cassert>
13 
14 template <
15   typename DerivedV,
16   typename DerivedEle,
17   typename Scalar>
normal_derivative(const Eigen::PlainObjectBase<DerivedV> & V,const Eigen::PlainObjectBase<DerivedEle> & Ele,Eigen::SparseMatrix<Scalar> & DD)18 IGL_INLINE void igl::normal_derivative(
19   const Eigen::PlainObjectBase<DerivedV> & V,
20   const Eigen::PlainObjectBase<DerivedEle> & Ele,
21   Eigen::SparseMatrix<Scalar>& DD)
22 {
23   using namespace Eigen;
24   using namespace std;
25   // Element simplex-size
26   const size_t ss = Ele.cols();
27   assert( ((ss==3) || (ss==4)) && "Only triangles or tets");
28   // cotangents
29   Matrix<Scalar,Dynamic,Dynamic> C;
30   cotmatrix_entries(V,Ele,C);
31   vector<Triplet<Scalar> > IJV;
32   // Number of elements
33   const size_t m = Ele.rows();
34   // Number of vertices
35   const size_t n = V.rows();
36   switch(ss)
37   {
38     default:
39       assert(false);
40       return;
41     case 4:
42     {
43       const MatrixXi DDJ =
44         slice(
45           Ele,
46           (VectorXi(24)<<
47             1,0,2,0,3,0,2,1,3,1,0,1,3,2,0,2,1,2,0,3,1,3,2,3).finished(),
48           2);
49       MatrixXi DDI(m,24);
50       for(size_t f = 0;f<4;f++)
51       {
52         const auto & I = (igl::LinSpaced<VectorXi >(m,0,m-1).array()+f*m).eval();
53         for(size_t r = 0;r<6;r++)
54         {
55           DDI.col(f*6+r) = I;
56         }
57       }
58       const DiagonalMatrix<Scalar,24,24> S =
59         (Matrix<Scalar,2,1>(1,-1).template replicate<12,1>()).asDiagonal();
60       Matrix<Scalar,Dynamic,Dynamic> DDV =
61         slice(
62           C,
63           (VectorXi(24)<<
64             2,2,1,1,3,3,0,0,4,4,2,2,5,5,1,1,0,0,3,3,4,4,5,5).finished(),
65           2);
66       DDV *= S;
67 
68       IJV.reserve(DDV.size());
69       for(size_t f = 0;f<6*4;f++)
70       {
71         for(size_t e = 0;e<m;e++)
72         {
73           IJV.push_back(Triplet<Scalar>(DDI(e,f),DDJ(e,f),DDV(e,f)));
74         }
75       }
76       DD.resize(m*4,n);
77       DD.setFromTriplets(IJV.begin(),IJV.end());
78       break;
79     }
80     case 3:
81     {
82       const MatrixXi DDJ =
83         slice(Ele,(VectorXi(12)<<2,0,1,0,0,1,2,1,1,2,0,2).finished(),2);
84       MatrixXi DDI(m,12);
85       for(size_t f = 0;f<3;f++)
86       {
87         const auto & I = (igl::LinSpaced<VectorXi >(m,0,m-1).array()+f*m).eval();
88         for(size_t r = 0;r<4;r++)
89         {
90           DDI.col(f*4+r) = I;
91         }
92       }
93       const DiagonalMatrix<Scalar,12,12> S =
94         (Matrix<Scalar,2,1>(1,-1).template replicate<6,1>()).asDiagonal();
95       Matrix<Scalar,Dynamic,Dynamic> DDV =
96         slice(C,(VectorXi(12)<<1,1,2,2,2,2,0,0,0,0,1,1).finished(),2);
97       DDV *= S;
98 
99       IJV.reserve(DDV.size());
100       for(size_t f = 0;f<12;f++)
101       {
102         for(size_t e = 0;e<m;e++)
103         {
104           IJV.push_back(Triplet<Scalar>(DDI(e,f),DDJ(e,f),DDV(e,f)));
105         }
106       }
107       DD.resize(m*3,n);
108       DD.setFromTriplets(IJV.begin(),IJV.end());
109       break;
110     }
111   }
112 
113 }
114 
115 #ifdef IGL_STATIC_LIBRARY
116 // Explicit template instantiation
117 template void igl::normal_derivative<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int>&);
118 #endif
119