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