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);
263 for (int i = 0; i < 4; i++)
264 {
265 const auto fid = M.elements[idx.element].fs[i];
266 const auto &fvid = M.faces[fid].vs;
267 int fv0 = fvid[0], fv1 = fvid[1], fv2 = fvid[2];
268 if (fv0 > fv2)
269 swap(fv0, fv2);
270 if (fv0 > fv1)
271 swap(fv0, fv1);
272 if (fv1 > fv2)
273 swap(fv1, fv2);
274
275 assert(fv0 < fv1);
276 assert(fv0 < fv2);
277 assert(fv1 < fv2);
278
279 if (v0 != fv0 || v1 != fv1 || v2 != fv2)
280 continue;
281
282 idx.face = fid;
283 idx.element_patch = i;
284
285 for (int j = 0; j < 3; j++)
286 {
287 const auto eid = M.faces[fid].es[j];
288 const auto &veid = M.edges[eid].vs;
289 assert(veid[0] < veid[1]);
290 if (veid[0] == v0_ && veid[1] == v1_)
291 {
292 idx.edge = eid;
293 if (fvid[0] == idx.vertex)
294 idx.face_corner = 0;
295 else if (fvid[1] == idx.vertex)
296 idx.face_corner = 1;
297 else
298 idx.face_corner = 2;
299
300 assert(idx.vertex == v0i);
301 assert(switch_vertex(M, idx).vertex == v1i);
302 assert(switch_vertex(M, switch_edge(M, idx)).vertex == v2i);
303 return idx;
304 }
305 }
306 }
307 }
308 assert(false);
309 return idx;
310}
311// Navigation in a surface mesh
312polyfem::mesh::Navigation3D::Index polyfem::mesh::Navigation3D::switch_vertex(const Mesh3DStorage &M, Index idx)
313{
314 // igl::Timer timer; timer.start();
315
316 if (M.type == MeshType::TET)
317 {
318 if (idx.vertex == M.EV(0, idx.edge))
319 idx.vertex = M.EV(1, idx.edge);
320 else
321 idx.vertex = M.EV(0, idx.edge);
322
323 if (M.FV(0, idx.face) == idx.vertex)
324 idx.face_corner = 0;
325 else if (M.FV(1, idx.face) == idx.vertex)
326 idx.face_corner = 1;
327 else
328 idx.face_corner = 2;
329 }
330 else
331 {
332 if (idx.vertex == M.edges[idx.edge].vs[0])
333 idx.vertex = M.edges[idx.edge].vs[1];
334 else
335 idx.vertex = M.edges[idx.edge].vs[0];
336
337 int &corner = idx.face_corner, n = M.faces[idx.face].vs.size(), corner_1 = (corner - 1 + n) % n, corner1 = (corner + 1) % n;
338 if (M.faces[idx.face].vs[corner1] == idx.vertex)
339 idx.face_corner = corner1;
340 else if (M.faces[idx.face].vs[corner_1] == idx.vertex)
341 idx.face_corner = corner_1;
342 }
343 // switch_vertex_time += timer.getElapsedTime();
344 return idx;
345}
346
347polyfem::mesh::Navigation3D::Index polyfem::mesh::Navigation3D::switch_edge(const Mesh3DStorage &M, Index idx)
348{
349 // igl::Timer timer; timer.start();
350
351 if (M.type == MeshType::TET)
352 {
353 if (idx.edge == M.FE(idx.face_corner, idx.face))
354 idx.edge = M.FE((idx.face_corner + 2) % 3, idx.face);
355 else
356 idx.edge = M.FE(idx.face_corner, idx.face);
357 }
358 else
359 {
360 int n = M.faces[idx.face].vs.size();
361 if (idx.edge == M.faces[idx.face].es[idx.face_corner])
362 idx.edge = M.faces[idx.face].es[(idx.face_corner - 1 + n) % n];
363 else
364 idx.edge = M.faces[idx.face].es[idx.face_corner];
365 }
366 // switch_edge_time += timer.getElapsedTime();
367 return idx;
368}
369
370polyfem::mesh::Navigation3D::Index polyfem::mesh::Navigation3D::switch_face(const Mesh3DStorage &M, Index idx)
371{
372 // igl::Timer timer; timer.start();
373 if (M.type == MeshType::TET)
374 {
375 for (int i = 0; i < 4; i++)
376 {
377 const auto &fid = M.HF(i, idx.element);
378 if (fid == idx.face)
379 continue;
380 if (M.FE(0, fid) == idx.edge || M.FE(1, fid) == idx.edge || M.FE(2, fid) == idx.edge)
381 {
382 idx.face = fid;
383 idx.element_patch = i;
384 if (M.FV(0, fid) == idx.vertex)
385 idx.face_corner = 0;
386 else if (M.FV(1, fid) == idx.vertex)
387 idx.face_corner = 1;
388 else
389 idx.face_corner = 2;
390 break;
391 }
392 }
393 }
394 else
395 {
396 const vector<uint32_t> &efs = M.edges[idx.edge].neighbor_fs, &hfs = M.elements[idx.element].fs;
397 std::array<uint32_t, 2> sharedfs;
398 int num = 2;
399 MeshProcessing3D::set_intersection_own(efs, hfs, sharedfs, num);
400 if (sharedfs[0] == idx.face)
401 idx.face = sharedfs[1];
402 else
403 idx.face = sharedfs[0];
404 for (int i = 0; i < hfs.size(); i++)
405 if (idx.face == hfs[i])
406 {
407 idx.element_patch = i;
408 break;
409 }
410
411 const vector<uint32_t> &fvs = M.faces[idx.face].vs;
412 for (int i = 0; i < fvs.size(); i++)
413 if (idx.vertex == fvs[i])
414 {
415 idx.face_corner = i;
416 break;
417 }
418 }
419
420 // switch_face_time += timer.getElapsedTime();
421
422 return idx;
423}
424
425polyfem::mesh::Navigation3D::Index polyfem::mesh::Navigation3D::switch_element(const Mesh3DStorage &M, Index idx)
426{
427 // igl::Timer timer; timer.start();
428
429 if (M.type == MeshType::TET)
430 {
431 if (M.FH(1, idx.face) == -1)
432 {
433 idx.element = -1;
434 return idx;
435 }
436 if (M.FH(0, idx.face) == idx.element)
437 {
438 idx.element = M.FH(1, idx.face);
439 idx.element_patch = M.FHi(1, idx.face);
440 }
441 else
442 {
443 idx.element = M.FH(0, idx.face);
444 idx.element_patch = M.FHi(0, idx.face);
445 }
446 }
447 else
448 {
449 if (M.faces[idx.face].neighbor_hs.size() == 1)
450 {
451 idx.element = -1;
452 return idx;
453 }
454 else
455 {
456 if (M.faces[idx.face].neighbor_hs[0] == idx.element)
457 idx.element = M.faces[idx.face].neighbor_hs[1];
458 else
459 idx.element = M.faces[idx.face].neighbor_hs[0];
460
461 const vector<uint32_t> &fs = M.elements[idx.element].fs;
462 for (int i = 0; i < fs.size(); i++)
463 if (idx.face == fs[i])
464 {
465 idx.element_patch = i;
466 break;
467 }
468 }
469 // const vector<uint32_t> &fvs = M.faces[idx.face].vs;
470 // for(int i=0;i<fvs.size();i++) if(idx.vertex == fvs[i]){idx.face_corner=i; break;}
471 }
472
473 // switch_element_time += timer.getElapsedTime();
474 return idx;
475}
polyfem::mesh::MeshProcessing3D::global_orientation_hexes
void global_orientation_hexes(Mesh3DStorage &hmi)
Definition MeshProcessing3D.cpp:679
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:1737
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:425
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:312
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:347
polyfem::mesh::Navigation3D::switch_face
Index switch_face(const Mesh3DStorage &M, Index idx)
Definition Navigation3D.cpp:370
polyfem::mesh::find
int find(const Eigen::VectorXi &vec, int x)
Definition NCMesh2D.cpp:289