PolyFEM
Loading...
Searching...
No Matches
Navigation3D.cpp
Go to the documentation of this file.
1#include "Navigation3D.hpp"
2#include "MeshProcessing3D.hpp"
3
4// #include <igl/Timer.h>
5#include <algorithm>
6#include <iterator>
7#include <set>
8#include <cassert>
9
10using namespace polyfem::mesh::Navigation3D;
11using namespace polyfem;
12using namespace std;
13
14// double polyfem::mesh::Navigation3D::get_index_from_element_face_time;
15// double polyfem::mesh::Navigation3D::switch_vertex_time;
16// double polyfem::mesh::Navigation3D::switch_edge_time;
17// double polyfem::mesh::Navigation3D::switch_face_time;
18// double polyfem::mesh::Navigation3D::switch_element_time;
19
20void polyfem::mesh::Navigation3D::prepare_mesh(Mesh3DStorage &M)
21{
22 if (M.type != MeshType::TET)
23 M.type = MeshType::HYB;
24 MeshProcessing3D::build_connectivity(M);
25 MeshProcessing3D::global_orientation_hexes(M);
26}
27
28polyfem::mesh::Navigation3D::Index polyfem::mesh::Navigation3D::get_index_from_element_face(const Mesh3DStorage &M, int hi)
29{
30 // igl::Timer timer; timer.start();
31
32 Index idx;
33 if (M.type == MeshType::TET)
34 {
35 idx.element = hi;
36 idx.element_patch = 0;
37 idx.face = M.HF(0, hi);
38 idx.face_corner = 0;
39 idx.vertex = M.FV(0, idx.face);
40 idx.edge = M.FE(0, idx.face);
41
42 if (M.elements[hi].fs_flag[idx.element_patch])
43 idx.edge = M.FE(2, idx.face);
44 // get_index_from_element_face_time += timer.getElapsedTime();
45 }
46 else if (M.elements[hi].hex)
47 {
48 idx.element = hi;
49 // idx.element_patch = 0;
50 // idx.face = M.elements[hi].fs[idx.element_patch];
51
52 // idx.vertex = M.elements[hi].vs[0];
53 // idx.face_corner = 0;
54 // idx.edge = M.faces[idx.face].es[0];
55
56 vector<uint32_t> fvs, fvs_;
57 fvs.insert(fvs.end(), M.elements[hi].vs.begin(), M.elements[hi].vs.begin() + 4);
58 sort(fvs.begin(), fvs.end());
59 idx.element_patch = -1;
60
61 for (uint32_t i = 0; i < 6; i++)
62 {
63 idx.element_patch = i;
64 fvs_ = M.faces[M.elements[hi].fs[i]].vs;
65 sort(fvs_.begin(), fvs_.end());
66 if (std::equal(fvs.begin(), fvs.end(), fvs_.begin()))
67 break;
68 }
69 idx.face = M.elements[hi].fs[idx.element_patch];
70
71 idx.vertex = M.elements[hi].vs[0];
72 idx.face_corner = find(M.faces[idx.face].vs.begin(), M.faces[idx.face].vs.end(), idx.vertex) - M.faces[idx.face].vs.begin();
73
74 int v0 = idx.vertex, v1 = M.elements[hi].vs[1];
75 const vector<uint32_t> &ves0 = M.vertices[v0].neighbor_es, &ves1 = M.vertices[v1].neighbor_es;
76 std::array<uint32_t, 2> sharedes;
77 int num = 1;
78 MeshProcessing3D::set_intersection_own(ves0, ves1, sharedes, num);
79 idx.edge = sharedes[0];
80 // get_index_from_element_face_time += timer.getElapsedTime();
81 }
82 else
83 {
84 idx = get_index_from_element_face(M, hi, 0, 0);
85 }
86
87 return idx;
88}
89
90polyfem::mesh::Navigation3D::Index polyfem::mesh::Navigation3D::get_index_from_element_face(const Mesh3DStorage &M, int hi, int lf, int lv)
91{
92 // igl::Timer timer; timer.start();
93 Index idx;
94
95 if (hi >= M.elements.size())
96 hi = hi % M.elements.size();
97 idx.element = hi;
98
99 if (lf >= M.elements[hi].fs.size())
100 lf = lf % M.elements[hi].fs.size();
101 idx.element_patch = lf;
102 idx.face = M.elements[hi].fs[idx.element_patch];
103
104 if (lv >= M.faces[idx.face].vs.size())
105 lv = lv % M.faces[idx.face].vs.size();
106 idx.face_corner = lv;
107 idx.vertex = M.faces[idx.face].vs[idx.face_corner];
108
109 int ei = idx.face_corner;
110 if (M.elements[hi].fs_flag[idx.element_patch])
111 ei = (idx.face_corner + M.faces[idx.face].vs.size() - 1) % M.faces[idx.face].vs.size();
112 idx.edge = M.faces[idx.face].es[ei];
113 // timer.stop();
114 // get_index_from_element_face_time += timer.getElapsedTime();
115
116 return idx;
117}
118polyfem::mesh::Navigation3D::Index polyfem::mesh::Navigation3D::get_index_from_element_edge(const Mesh3DStorage &M, int hi, int v0i, int v1i)
119{
120 Index idx;
121 idx.element = hi;
122 idx.vertex = v0i;
123 int v0 = v0i;
124 int v1 = v1i;
125 if (v0 > v1)
126 std::swap(v0, v1);
127 assert(v0 < v1);
128
129 if (M.type == MeshType::TET)
130 {
131 for (int i = 0; i < 4; i++)
132 {
133 const auto &fid = M.HF(i, idx.element);
134 for (int j = 0; j < 3; j++)
135 {
136 const auto &eid = M.FE(j, fid);
137 assert(M.EV(0, eid) < M.EV(1, eid));
138 if (M.EV(0, eid) == v0 && M.EV(1, eid) == v1)
139 {
140 idx.element_patch = i;
141 idx.face = fid;
142 idx.edge = eid;
143 if (M.FV(0, fid) == idx.vertex)
144 idx.face_corner = 0;
145 else if (M.FV(1, fid) == idx.vertex)
146 idx.face_corner = 1;
147 else
148 idx.face_corner = 2;
149
150 assert(idx.vertex == v0i);
151 assert(switch_vertex(M, idx).vertex == v1i);
152 assert(idx.element == hi);
153
154 return idx;
155 }
156 }
157 }
158 }
159 else
160 {
161 for (int i = 0; i < M.elements[hi].fs.size(); i++)
162 {
163 const auto &fid = M.elements[hi].fs[i];
164 for (int j = 0; j < M.faces[fid].es.size(); j++)
165 {
166 const auto &eid = M.faces[fid].es[j];
167 assert(M.edges[eid].vs[0] < M.edges[eid].vs[1]);
168 if (M.edges[eid].vs[0] == v0 && M.edges[eid].vs[1] == v1)
169 {
170 idx.element_patch = i;
171 idx.face = fid;
172 idx.edge = eid;
173 for (int k = 0; k < M.faces[fid].vs.size(); k++)
174 if (M.faces[fid].vs[k] == idx.vertex)
175 idx.face_corner = k;
176
177 assert(idx.vertex == v0i);
178 assert(switch_vertex(M, idx).vertex == v1i);
179 assert(idx.element == hi);
180
181 return idx;
182 }
183 }
184 }
185 }
186
187 assert(false);
188 return idx;
189}
190
191polyfem::mesh::Navigation3D::Index polyfem::mesh::Navigation3D::get_index_from_element_tri(const Mesh3DStorage &M, int hi, int v0i, int v1i, int v2i)
192{
193 int v0 = v0i;
194 int v1 = v1i;
195 int v2 = v2i;
196
197 Index idx;
198 idx.element = hi;
199 idx.vertex = v0;
200 int v0_ = v0, v1_ = v1;
201 if (v0_ > v1_)
202 swap(v0_, v1_);
203
204 if (v0 > v2)
205 swap(v0, v2);
206 if (v0 > v1)
207 swap(v0, v1);
208 if (v1 > v2)
209 swap(v1, v2);
210
211 assert(v0 < v1);
212 assert(v0 < v2);
213 assert(v1 < v2);
214
215 if (M.type == MeshType::TET)
216 {
217 for (int i = 0; i < 4; i++)
218 {
219 const auto &fid = M.HF(i, idx.element);
220 int fv0 = M.FV(0, fid), fv1 = M.FV(1, fid), fv2 = M.FV(2, fid);
221 if (fv0 > fv2)
222 swap(fv0, fv2);
223 if (fv0 > fv1)
224 swap(fv0, fv1);
225 if (fv1 > fv2)
226 swap(fv1, fv2);
227
228 assert(fv0 < fv1);
229 assert(fv0 < fv2);
230 assert(fv1 < fv2);
231
232 if (v0 != fv0 || v1 != fv1 || v2 != fv2)
233 continue;
234
235 idx.face = fid;
236 idx.element_patch = i;
237
238 for (int j = 0; j < 3; j++)
239 {
240 const auto &eid = M.FE(j, fid);
241 assert(M.EV(0, eid) < M.EV(1, eid));
242 if (M.EV(0, eid) == v0_ && M.EV(1, eid) == v1_)
243 {
244 idx.edge = eid;
245 if (M.FV(0, fid) == idx.vertex)
246 idx.face_corner = 0;
247 else if (M.FV(1, fid) == idx.vertex)
248 idx.face_corner = 1;
249 else
250 idx.face_corner = 2;
251
252 assert(idx.vertex == v0i);
253 assert(switch_vertex(M, idx).vertex == v1i);
254 assert(switch_vertex(M, switch_edge(M, idx)).vertex == v2i);
255 return idx;
256 }
257 }
258 }
259 }
260 else
261 {
262 assert(M.elements[idx.element].fs.size() == 4 || M.elements[idx.element].fs.size() == 5);
263 for (int i = 0; i < M.elements[idx.element].fs.size(); i++)
264 {
265 const auto fid = M.elements[idx.element].fs[i];
266 const auto &fvid = M.faces[fid].vs;
267 if (fvid.size() != 3)
268 continue;
269 int fv0 = fvid[0], fv1 = fvid[1], fv2 = fvid[2];
270 if (fv0 > fv2)
271 swap(fv0, fv2);
272 if (fv0 > fv1)
273 swap(fv0, fv1);
274 if (fv1 > fv2)
275 swap(fv1, fv2);
276
277 assert(fv0 < fv1);
278 assert(fv0 < fv2);
279 assert(fv1 < fv2);
280
281 if (v0 != fv0 || v1 != fv1 || v2 != fv2)
282 continue;
283
284 idx.face = fid;
285 idx.element_patch = i;
286
287 for (int j = 0; j < 3; j++)
288 {
289 const auto eid = M.faces[fid].es[j];
290 const auto &veid = M.edges[eid].vs;
291 assert(veid[0] < veid[1]);
292 if (veid[0] == v0_ && veid[1] == v1_)
293 {
294 idx.edge = eid;
295 if (fvid[0] == idx.vertex)
296 idx.face_corner = 0;
297 else if (fvid[1] == idx.vertex)
298 idx.face_corner = 1;
299 else
300 idx.face_corner = 2;
301
302 assert(idx.vertex == v0i);
303 assert(switch_vertex(M, idx).vertex == v1i);
304 assert(switch_vertex(M, switch_edge(M, idx)).vertex == v2i);
305 return idx;
306 }
307 }
308 }
309 }
310 assert(false);
311 return idx;
312}
313// Navigation in a surface mesh
314polyfem::mesh::Navigation3D::Index polyfem::mesh::Navigation3D::switch_vertex(const Mesh3DStorage &M, Index idx)
315{
316 // igl::Timer timer; timer.start();
317
318 if (M.type == MeshType::TET)
319 {
320 if (idx.vertex == M.EV(0, idx.edge))
321 idx.vertex = M.EV(1, idx.edge);
322 else
323 idx.vertex = M.EV(0, idx.edge);
324
325 if (M.FV(0, idx.face) == idx.vertex)
326 idx.face_corner = 0;
327 else if (M.FV(1, idx.face) == idx.vertex)
328 idx.face_corner = 1;
329 else
330 idx.face_corner = 2;
331 }
332 else
333 {
334 if (idx.vertex == M.edges[idx.edge].vs[0])
335 idx.vertex = M.edges[idx.edge].vs[1];
336 else
337 idx.vertex = M.edges[idx.edge].vs[0];
338
339 int &corner = idx.face_corner, n = M.faces[idx.face].vs.size(), corner_1 = (corner - 1 + n) % n, corner1 = (corner + 1) % n;
340 if (M.faces[idx.face].vs[corner1] == idx.vertex)
341 idx.face_corner = corner1;
342 else if (M.faces[idx.face].vs[corner_1] == idx.vertex)
343 idx.face_corner = corner_1;
344 }
345 // switch_vertex_time += timer.getElapsedTime();
346 return idx;
347}
348
349polyfem::mesh::Navigation3D::Index polyfem::mesh::Navigation3D::switch_edge(const Mesh3DStorage &M, Index idx)
350{
351 // igl::Timer timer; timer.start();
352
353 if (M.type == MeshType::TET)
354 {
355 if (idx.edge == M.FE(idx.face_corner, idx.face))
356 idx.edge = M.FE((idx.face_corner + 2) % 3, idx.face);
357 else
358 idx.edge = M.FE(idx.face_corner, idx.face);
359 }
360 else
361 {
362 int n = M.faces[idx.face].vs.size();
363 if (idx.edge == M.faces[idx.face].es[idx.face_corner])
364 idx.edge = M.faces[idx.face].es[(idx.face_corner - 1 + n) % n];
365 else
366 idx.edge = M.faces[idx.face].es[idx.face_corner];
367 }
368 // switch_edge_time += timer.getElapsedTime();
369 return idx;
370}
371
372polyfem::mesh::Navigation3D::Index polyfem::mesh::Navigation3D::switch_face(const Mesh3DStorage &M, Index idx)
373{
374 // igl::Timer timer; timer.start();
375 if (M.type == MeshType::TET)
376 {
377 for (int i = 0; i < 4; i++)
378 {
379 const auto &fid = M.HF(i, idx.element);
380 if (fid == idx.face)
381 continue;
382 if (M.FE(0, fid) == idx.edge || M.FE(1, fid) == idx.edge || M.FE(2, fid) == idx.edge)
383 {
384 idx.face = fid;
385 idx.element_patch = i;
386 if (M.FV(0, fid) == idx.vertex)
387 idx.face_corner = 0;
388 else if (M.FV(1, fid) == idx.vertex)
389 idx.face_corner = 1;
390 else
391 idx.face_corner = 2;
392 break;
393 }
394 }
395 }
396 else
397 {
398 const vector<uint32_t> &efs = M.edges[idx.edge].neighbor_fs, &hfs = M.elements[idx.element].fs;
399 std::array<uint32_t, 2> sharedfs;
400 int num = 2;
401 MeshProcessing3D::set_intersection_own(efs, hfs, sharedfs, num);
402 if (sharedfs[0] == idx.face)
403 idx.face = sharedfs[1];
404 else
405 idx.face = sharedfs[0];
406 for (int i = 0; i < hfs.size(); i++)
407 if (idx.face == hfs[i])
408 {
409 idx.element_patch = i;
410 break;
411 }
412
413 const vector<uint32_t> &fvs = M.faces[idx.face].vs;
414 for (int i = 0; i < fvs.size(); i++)
415 if (idx.vertex == fvs[i])
416 {
417 idx.face_corner = i;
418 break;
419 }
420 }
421
422 // switch_face_time += timer.getElapsedTime();
423
424 return idx;
425}
426
427polyfem::mesh::Navigation3D::Index polyfem::mesh::Navigation3D::switch_element(const Mesh3DStorage &M, Index idx)
428{
429 // igl::Timer timer; timer.start();
430
431 if (M.type == MeshType::TET)
432 {
433 if (M.FH(1, idx.face) == -1)
434 {
435 idx.element = -1;
436 return idx;
437 }
438 if (M.FH(0, idx.face) == idx.element)
439 {
440 idx.element = M.FH(1, idx.face);
441 idx.element_patch = M.FHi(1, idx.face);
442 }
443 else
444 {
445 idx.element = M.FH(0, idx.face);
446 idx.element_patch = M.FHi(0, idx.face);
447 }
448 }
449 else
450 {
451 if (M.faces[idx.face].neighbor_hs.size() == 1)
452 {
453 idx.element = -1;
454 return idx;
455 }
456 else
457 {
458 if (M.faces[idx.face].neighbor_hs[0] == idx.element)
459 idx.element = M.faces[idx.face].neighbor_hs[1];
460 else
461 idx.element = M.faces[idx.face].neighbor_hs[0];
462
463 const vector<uint32_t> &fs = M.elements[idx.element].fs;
464 for (int i = 0; i < fs.size(); i++)
465 if (idx.face == fs[i])
466 {
467 idx.element_patch = i;
468 break;
469 }
470 }
471 // const vector<uint32_t> &fvs = M.faces[idx.face].vs;
472 // for(int i=0;i<fvs.size();i++) if(idx.vertex == fvs[i]){idx.face_corner=i; break;}
473 }
474
475 // switch_element_time += timer.getElapsedTime();
476 return idx;
477}
polyfem::mesh::MeshProcessing3D::global_orientation_hexes
void global_orientation_hexes(Mesh3DStorage &hmi)
Definition MeshProcessing3D.cpp:675
polyfem::mesh::MeshProcessing3D::set_intersection_own
void set_intersection_own(const std::vector< uint32_t > &A, const std::vector< uint32_t > &B, std::array< uint32_t, 2 > &C, int &num)
Definition MeshProcessing3D.cpp:1718
polyfem::mesh::MeshProcessing3D::build_connectivity
void build_connectivity(Mesh3DStorage &hmi)
Definition MeshProcessing3D.cpp:18
polyfem::mesh::Navigation3D::switch_element
Index switch_element(const Mesh3DStorage &M, Index idx)
Definition Navigation3D.cpp:427
polyfem::mesh::Navigation3D::get_index_from_element_tri
Index get_index_from_element_tri(const Mesh3DStorage &M, int hi, int v0, int v1, int v2)
Definition Navigation3D.cpp:191
polyfem::mesh::Navigation3D::switch_vertex
Index switch_vertex(const Mesh3DStorage &M, Index idx)
Definition Navigation3D.cpp:314
polyfem::mesh::Navigation3D::get_index_from_element_face
Index get_index_from_element_face(const Mesh3DStorage &M, int hi)
Definition Navigation3D.cpp:28
polyfem::mesh::Navigation3D::get_index_from_element_edge
Index get_index_from_element_edge(const Mesh3DStorage &M, int hi, int v0, int v1)
Definition Navigation3D.cpp:118
polyfem::mesh::Navigation3D::prepare_mesh
void prepare_mesh(Mesh3DStorage &M)
Definition Navigation3D.cpp:20
polyfem::mesh::Navigation3D::switch_edge
Index switch_edge(const Mesh3DStorage &M, Index idx)
Definition Navigation3D.cpp:349
polyfem::mesh::Navigation3D::switch_face
Index switch_face(const Mesh3DStorage &M, Index idx)
Definition Navigation3D.cpp:372
polyfem::mesh::find
int find(const Eigen::VectorXi &vec, int x)
Definition NCMesh2D.cpp:289