10 const int dim = mesh.dimension();
11 const int n_verts = mesh.n_vertices();
12 assert(mesh.is_simplicial());
14 std::vector<int> tmp = {};
15 std::set<int> surface_ids = std::set(tmp.begin(), tmp.end());
18 std::vector<bool> active_mask;
19 active_mask.assign(n_verts,
false);
20 std::vector<Eigen::Triplet<bool>> T_adj;
22 for (
int b = 0; b < mesh.n_boundary_elements(); b++)
24 const int boundary_id = mesh.get_boundary_id(b);
25 if (!surface_ids.empty() && surface_ids.find(boundary_id) == surface_ids.end())
28 for (
int lv = 0; lv < dim; lv++)
30 active_mask[mesh.boundary_element_vertex(b, lv)] =
true;
33 for (
int lv1 = 0; lv1 < dim; lv1++)
34 for (
int lv2 = 0; lv2 < lv1; lv2++)
36 T_adj.emplace_back(mesh.boundary_element_vertex(b, lv2), mesh.boundary_element_vertex(b, lv1),
true);
37 T_adj.emplace_back(mesh.boundary_element_vertex(b, lv1), mesh.boundary_element_vertex(b, lv2),
true);
42 adj.resize(n_verts, n_verts);
43 adj.setFromTriplets(T_adj.begin(), T_adj.end());
45 std::vector<int> degrees(n_verts, 0);
46 for (
int k = 0; k <
adj.outerSize(); ++k)
47 for (Eigen::SparseMatrix<bool, Eigen::RowMajor>::InnerIterator it(
adj, k); it; ++it)
51 L.resize(n_verts, n_verts);
54 std::vector<Eigen::Triplet<double>> T_L;
55 for (
int k = 0; k <
adj.outerSize(); ++k)
59 T_L.emplace_back(k, k, degrees[k]);
60 for (Eigen::SparseMatrix<bool, Eigen::RowMajor>::InnerIterator it(
adj, k); it; ++it)
62 assert(it.row() == k);
63 T_L.emplace_back(it.row(), it.col(), -1);
66 L.setFromTriplets(T_L.begin(), T_L.end());
67 L.prune([](
int i,
int j,
double val) {
return abs(
val) > 1e-12; });
115 const int dim = mesh.dimension();
116 const int n_verts = mesh.n_vertices();
118 Eigen::VectorXd grad;
121 grad.setZero(n_verts * dim);
122 for (
int b = 0; b <
adj.rows(); b++)
127 auto sum_normalized = s;
129 for (Eigen::SparseMatrix<bool, Eigen::RowMajor>::InnerIterator it(
adj, b); it; ++it)
131 assert(it.col() != b);
132 auto x = mesh.point(b) - mesh.point(it.col());
134 sum_norm +=
x.norm();
135 sum_normalized +=
x.normalized();
141 const double coeff =
power_ * pow(s.norm(),
power_ - 2.) / sum_norm;
143 grad.segment(b * dim, dim) += (s * valence - s.squaredNorm() * sum_normalized) * coeff;
144 for (Eigen::SparseMatrix<bool, Eigen::RowMajor>::InnerIterator it(
adj, b); it; ++it)
145 grad.segment(it.col() * dim, dim) -= (s + s.squaredNorm() * (mesh.point(it.col()) - mesh.point(b)).normalized()) * coeff;