generated from libigl/libigl-example-project
-
Notifications
You must be signed in to change notification settings - Fork 0
/
cotmatrix_timed.cpp
105 lines (94 loc) · 4.42 KB
/
cotmatrix_timed.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
// 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 "cotmatrix_timed.h"
#include <vector>
// For error printing
#include <cstdio>
#include <igl/cotmatrix_entries.h>
// Bug in unsupported/Eigen/SparseExtra needs iostream first
#include <iostream>
#include <Eigen/Dense>
#include <Eigen/Sparse>
template <typename DerivedV, typename DerivedF, typename Scalar>
IGL_INLINE void cotmatrix_timed(
const Eigen::MatrixBase<DerivedV> & V,
const Eigen::MatrixBase<DerivedF> & F,
Eigen::SparseMatrix<Scalar>& L)
{
using namespace Eigen;
using namespace std;
L.resize(V.rows(),V.rows());
Matrix<int,Dynamic,2> 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
L.reserve(10*V.rows());
edges.resize(3,2);
edges <<
1,2,
2,0,
0,1;
}else if(simplex_size == 4)
{
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<Scalar,Dynamic,Dynamic> C;
igl::cotmatrix_entries(V,F,C);
vector<Triplet<Scalar> > IJV;
chrono::steady_clock::time_point begin_pushback = chrono::steady_clock::now();
IJV.reserve(F.rows()*edges.rows()*4);
// Loop over triangles
for(int i = 0; i < F.rows(); i++)
{
// loop over edges of element
for(int e = 0;e<edges.rows();e++)
{
int source = F(i,edges(e,0));
int dest = F(i,edges(e,1));
IJV.emplace_back(source,dest,C(i,e));
IJV.emplace_back(dest,source,C(i,e));
IJV.emplace_back(source,source,-C(i,e));
IJV.emplace_back(dest,dest,-C(i,e));
}
}
chrono::steady_clock::time_point end_pushback = chrono::steady_clock::now();
std::chrono::duration<double, std::milli> time_pushback = end_pushback-begin_pushback;
chrono::steady_clock::time_point begin_set = chrono::steady_clock::now();
L.setFromTriplets(IJV.begin(),IJV.end());
chrono::steady_clock::time_point end_set = chrono::steady_clock::now();
std::chrono::duration<double, std::milli> time_set = end_set-begin_set;
std::chrono::duration<double, std::milli> time_assembly = end_set-begin_pushback;
cout << time_assembly.count() << " " << std::flush;
}
#ifdef IGL_STATIC_LIBRARY
// Explicit template instantiation
// generated by autoexplicit.sh
template void cotmatrix_timed<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -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::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::SparseMatrix<double, 0, int>&, Eigen::SparseMatrix<double, 0, int>&, Eigen::SparseMatrix<double, 0, int>&);
// generated by autoexplicit.sh
template void cotmatrix_timed<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&, Eigen::SparseMatrix<double, 0, int>&);
template void cotmatrix_timed<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&, Eigen::SparseMatrix<double, 0, int>&);
template void cotmatrix_timed<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&, Eigen::SparseMatrix<double, 0, int>&);
template void cotmatrix_timed<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>&);
#endif