PolyFEM
Loading...
Searching...
No Matches
Singularities.cpp
Go to the documentation of this file.
1
2#include "Singularities.hpp"
3#include "Navigation.hpp"
5#include <algorithm>
7
9 const GEO::Mesh &M, Eigen::VectorXi &V, int regular_degree, bool ignore_border)
10{
11 using GEO::index_t;
12 std::vector<int> degree(M.vertices.nb(), 0);
13 for (index_t f = 0; f < M.facets.nb(); ++f)
14 {
15 for (index_t lv = 0; lv < M.facets.nb_vertices(f); ++lv)
16 {
17 index_t v = M.facets.vertex(f, lv);
18 degree[v]++;
19 }
20 }
21
22 // Ignore border vertices if requested
23 if (ignore_border)
24 {
25 GEO::Attribute<bool> boundary_vertices(M.vertices.attributes(), "boundary_vertex");
26 for (index_t v = 0; v < M.vertices.nb(); ++v)
27 {
28 if (boundary_vertices[v])
29 {
30 degree[v] = regular_degree;
31 }
32 }
33 }
34
35 // Count and create list of irregular vertices
36 long nb_reg = std::count(degree.begin(), degree.end(), regular_degree);
37 V.resize(M.vertices.nb() - nb_reg);
38 for (index_t v = 0, i = 0; v < M.vertices.nb(); ++v)
39 {
40 if (degree[v] != regular_degree)
41 {
42 assert(i < V.size());
43 V(i++) = v;
44 }
45 }
46}
47
48void polyfem::mesh::singular_edges(const GEO::Mesh &M, const Eigen::VectorXi &V, Eigen::MatrixX2i &E)
49{
50 using GEO::index_t;
51 std::vector<bool> is_singular(M.vertices.nb(), false);
52 for (int i = 0; i < V.size(); ++i)
53 {
54 is_singular[V(i)] = true;
55 }
56 std::vector<std::pair<int, int>> edges;
57 for (index_t e = 0; e < M.edges.nb(); ++e)
58 {
59 index_t v1 = M.edges.vertex(e, 0);
60 index_t v2 = M.edges.vertex(e, 1);
61 if (is_singular[v1] && is_singular[v2])
62 {
63 edges.emplace_back(v1, v2);
64 }
65 }
66 E.resize(edges.size(), 2);
67 int cnt = 0;
68 for (const auto &kv : edges)
69 {
70 assert(cnt < E.rows());
71 E.row(cnt++) << kv.first, kv.second;
72 }
73}
74
76 const GEO::Mesh &M, Eigen::VectorXi &V, Eigen::MatrixX2i &E, int regular_degree, bool ignore_border)
77{
78 singular_vertices(M, V, regular_degree, ignore_border);
79 singular_edges(M, V, E);
80}
81
83
85 GEO::Mesh &M, const Eigen::VectorXi &V, const Eigen::MatrixX2i &E, double t)
86{
87 using GEO::index_t;
89
90 assert(M.vertices.dimension() == 2 || M.vertices.dimension() == 3);
91
92 // Step 0: Make sure all endpoints of E are in V
93 std::vector<int> singulars;
94 singulars.reserve(V.size() + E.size());
95 for (int i = 0; i < V.size(); ++i)
96 {
97 singulars.push_back(V(i));
98 }
99 for (int i = 0; i < E.size(); ++i)
100 {
101 singulars.push_back(E.data()[i]);
102 }
103 std::sort(singulars.begin(), singulars.end());
104 auto it = std::unique(singulars.begin(), singulars.end());
105 singulars.resize(std::distance(singulars.begin(), it));
106
107 std::vector<bool> is_singular(M.vertices.nb(), false);
108 for (int v : singulars)
109 {
110 is_singular[v] = true;
111 }
112
113 // Step 1: Find all edges around singularities and create new vertices on those edge
114 std::vector<int> edge_to_midpoint(M.edges.nb(), -1);
115 for (index_t e = 0; e < M.edges.nb(); ++e)
116 {
117 for (index_t lv = 0; lv < 2; ++lv)
118 {
119 int v1 = M.edges.vertex(e, lv);
120 int v2 = M.edges.vertex(e, (lv + 1) % 2);
121 if (is_singular[v1] && !is_singular[v2])
122 {
123 GEO::vec3 coords = t * mesh_vertex(M, v1) + (1.0 - t) * mesh_vertex(M, v2);
124 edge_to_midpoint[e] = mesh_create_vertex(M, coords);
125 }
126 }
127 }
128
129 GEO::Attribute<GEO::index_t> c2e(M.facet_corners.attributes(), "edge_id");
130
131 // Step 2: Iterate over all facets, and subdivide accordingly
132 std::vector<index_t> facets_to_delete;
133 std::vector<int> next_vertex_around_hole(M.vertices.nb(), -1);
134 for (index_t f = 0, old_num_facets = M.facets.nb(); f < old_num_facets; ++f)
135 {
136 // a. Iterate over edges around the facet, until one marked edge is found
137 // b. Place the navigator index so that its vertex points to the unmarked vertex
138 // c. Navigate around until another marked edge is found (assert it is different from the previous one)
139 // d. Possibly invert the list of vertices / edge midpoints to respect initial facet orientation
140 // e. Create new facets given the list + 2 edge midpoints. The rules are:
141 // i. If #v == 1, don't create the midfacet vertex iff the original facet was a triangle
142 // ii. If #v == 2, then don't create the midfacet vertex
143 // iii. If #v >= 3, then create midfacet vertex iff orignal facet was a quad or a tri
144 assert(M.facets.nb_vertices(f) > 2);
145
146 // a. Iterate over edges around the facet, until one marked edge is found
147 Index idx1 = Navigation::get_index_from_face(M, c2e, f, 0);
148
149 for (index_t lv = 0; lv < M.facets.nb_vertices(f); ++lv)
150 {
151 if (edge_to_midpoint[idx1.edge] != -1)
152 {
153 break;
154 }
155 idx1 = Navigation::next_around_face(M, c2e, idx1);
156 }
157 if (edge_to_midpoint[idx1.edge] == -1)
158 {
159 continue;
160 }
161 facets_to_delete.push_back(f);
162
163 // b. Place the navigator index so that its vertex points to the unmarked vertex
164 if (is_singular[idx1.vertex])
165 {
166 idx1 = Navigation::switch_vertex(M, idx1);
167 assert(!is_singular[idx1.vertex]);
168 }
169
170 // c. Navigate around until another marked edge is found (assert it is different from the previous one)
171 Index idx2 = Navigation::switch_edge(M, c2e, idx1);
172 GEO::vector<index_t> unmarked_vertices;
173 for (index_t lv = 0; lv < M.facets.nb_vertices(f); ++lv)
174 {
175 assert(!is_singular[idx2.vertex]);
176 unmarked_vertices.push_back(idx2.vertex);
177 if (edge_to_midpoint[idx2.edge] != -1)
178 {
179 break;
180 }
181 idx2 = Navigation::next_around_face(M, c2e, idx2);
182 }
183 assert(edge_to_midpoint[idx2.edge] != -1);
184 assert(idx1.edge != idx2.edge);
185 assert(unmarked_vertices.front() == idx1.vertex);
186 assert(unmarked_vertices.back() == idx2.vertex);
187
188 // d. Possibly invert the list of vertices / edge midpoints to respect initial facet orientation
189 {
190 int v1 = idx2.vertex;
191 int v0 = Navigation::switch_vertex(M, idx2).vertex;
192 int lv1 = M.facets.find_vertex(f, v1);
193 int v2 = M.facets.vertex(f, (lv1 + 1) % M.facets.nb_vertices(f));
194 if (v0 != v2)
195 {
196 // Need to reverse the sequence in `unmarked_vertices`, and swap(idx1, idx2)
197 std::reverse(unmarked_vertices.begin(), unmarked_vertices.end());
198 std::swap(idx1, idx2);
199 }
200 }
201
202 // e. Create new facets given the list + 2 edge midpoints. The rules are:
203 // i. If #v == 1, don't create the midfacet vertex iff the original facet was a triangle
204 // ii. If #v == 2, then don't create the midfacet vertex
205 // iii. If #v >= 3, then create midfacet vertex iff orignal facet was a quad or a tri
206 int sf = M.facets.nb_vertices(f);
207 int nv = (int)unmarked_vertices.size();
208 bool create_midfacet_point = false;
209 if (nv == 1)
210 {
211 if (sf == 3)
212 {
213 create_midfacet_point = false;
214 }
215 else
216 {
217 create_midfacet_point = true;
218 }
219 }
220 else if (nv == 2)
221 {
222 create_midfacet_point = false;
223 }
224 else
225 {
226 if (sf == 3 || sf == 4)
227 {
228 create_midfacet_point = true;
229 }
230 else
231 {
232 create_midfacet_point = false;
233 }
234 }
235
236 int ve1 = edge_to_midpoint[idx1.edge];
237 int ve2 = edge_to_midpoint[idx2.edge];
238 if (create_midfacet_point)
239 {
240 GEO::vec3 coords = facet_barycenter(M, f);
241 index_t vf = mesh_create_vertex(M, coords);
242 next_vertex_around_hole.resize(vf + 1);
243 assert(next_vertex_around_hole[ve1] == -1);
244 next_vertex_around_hole[ve1] = vf;
245 next_vertex_around_hole[vf] = ve2;
246 if (nv == 1)
247 {
248 M.facets.create_quad(unmarked_vertices.front(), ve2, vf, ve1);
249 }
250 else
251 {
252 M.facets.create_quad(unmarked_vertices[0], unmarked_vertices[1], vf, ve1);
253 M.facets.create_quad(*(unmarked_vertices.rbegin() + 1), unmarked_vertices.back(), ve2, vf);
254 for (int i = 1; i < int(unmarked_vertices.size()) - 2; ++i)
255 {
256 M.facets.create_triangle(unmarked_vertices[i], unmarked_vertices[i + 1], vf);
257 }
258 }
259 }
260 else
261 {
262 assert(next_vertex_around_hole[ve1] == -1);
263 next_vertex_around_hole[ve1] = ve2;
264 unmarked_vertices.push_back(ve2);
265 unmarked_vertices.push_back(ve1);
266 M.facets.create_polygon(unmarked_vertices);
267 }
268 }
269
270 // Step 3: Iterate over new vertices, loop around holes and create polygonal facets
271 std::vector<bool> marked(M.vertices.nb(), false);
272 for (index_t v = 0; v < M.vertices.nb(); ++v)
273 {
274 if (next_vertex_around_hole[v] == -1 || marked[v])
275 {
276 continue;
277 }
278 GEO::vector<index_t> hole;
279 int vi = v;
280 do
281 {
282 vi = next_vertex_around_hole[vi];
283 assert(vi != -1);
284 hole.push_back(vi);
285 marked[vi] = true;
286 } while (vi != (int)v);
287 M.facets.create_polygon(hole);
288 }
289
290 // Step 4: Remove facets that have been split, and delete isolated vertices
291 GEO::vector<index_t> to_delete_mask(M.facets.nb(), 0u);
292 for (index_t f : facets_to_delete)
293 {
294 to_delete_mask[f] = 1;
295 }
296 M.facets.delete_elements(to_delete_mask);
297 M.edges.clear();
298 M.vertices.remove_isolated();
299
300#if 0
301 // Step 2: Find all faces adjacent to a marked edge and place a new vertex at the center of the facet if needed
302 std::vector<int> facet_to_midpoint(M.facets.nb(), -1);
303 GEO::vector<index_t> facets_to_delete;
304 for (index_t f = 0; f < M.facets.nb(); ++f) {
305 Index idx = Navigation::get_index_from_face(M, f, 0);
306 std::vector<int> marked_edges;
307 for (index_t le = 0; le < M.facets.nb_vertices(f); ++le) {
308 if (edge_to_midpoint[idx.edge] != -1) {
309 marked_edges.push_back(le);
310 }
311 idx = Navigation::next_around_face(M, idx);
312 }
313 // There should be exactly 2 split edges per face, otherwise I don't know what to do!
314 assert(marked_edges.empty() || marked_edges.size() == 2);
315 if (!marked_edges.empty()) {
316 int e1 = marked_edges.front();
317 int e2 = marked_edges.back();
318 if (e1 + 1 == e2 || (e2 + 1) % (int) M.facets.nb_vertices(f) == e1) {
319 GEO::vec3 coords = facet_barycenter(f);
320 facet_to_midpoint[f] = M.vertices.create_vertex(&coords[0]);
321 }
322 facets_to_delete.push_back(f);
323 }
324 }
325
326 // Step 3: Create new facets for each marked facet, based on its configuration
327 std::vector<index_t> next_around_hole(M.vertices.nb(), -1);
328 std::vector<index_t> prev_around_hole(M.vertices.nb(), -1);
329
330 // Step 4: Remove facets that have been split, and delete isolated vertices
331 M.facets.delete_elements(facets_to_delete);
332 M.vertices.remove_isolated();
333 // + orient facets
334#endif
335}
int V
void create_patch_around_singularities(GEO::Mesh &M, const Eigen::VectorXi &V, const Eigen::MatrixX2i &E, double t=0.5)
GEO::vec3 mesh_vertex(const GEO::Mesh &M, GEO::index_t v)
Retrieve a 3D vector with the position of a given vertex.
Definition MeshUtils.cpp:44
GEO::index_t mesh_create_vertex(GEO::Mesh &M, const GEO::vec3 &p)
Definition MeshUtils.cpp:77
void singularity_graph(const GEO::Mesh &M, Eigen::VectorXi &V, Eigen::MatrixX2i &E, int regular_degree=4, bool ignore_border=true)
GEO::vec3 facet_barycenter(const GEO::Mesh &M, GEO::index_t f)
Definition MeshUtils.cpp:64
void singular_edges(const GEO::Mesh &M, const Eigen::VectorXi &V, Eigen::MatrixX2i &E)
void singular_vertices(const GEO::Mesh &M, Eigen::VectorXi &V, int regular_degree=4, bool ignore_border=true)