PolyFEM
Loading...
Searching...
No Matches
MeshProcessing3D.cpp
Go to the documentation of this file.
1#include "MeshProcessing3D.hpp"
2#include <polyfem/utils/Logger.hpp>
3
4#include <Eigen/Dense>
5
6#include <algorithm>
7#include <map>
8#include <set>
9#include <queue>
10#include <iterator>
11#include <cassert>
12#include <limits>
13
14using namespace polyfem::mesh;
15using namespace polyfem;
16using namespace std;
17using namespace Eigen;
18
19void MeshProcessing3D::build_connectivity(Mesh3DStorage &hmi)
20{
21 hmi.edges.clear();
22 if (hmi.type == MeshType::TRI || hmi.type == MeshType::QUA || hmi.type == MeshType::H_SUR)
23 {
24 std::vector<std::tuple<uint32_t, uint32_t, uint32_t, uint32_t>> temp;
25 temp.reserve(hmi.faces.size() * 3);
26 for (uint32_t i = 0; i < hmi.faces.size(); ++i)
27 {
28 int vn = hmi.faces[i].vs.size();
29 for (uint32_t j = 0; j < vn; ++j)
30 {
31 uint32_t v0 = hmi.faces[i].vs[j], v1 = hmi.faces[i].vs[(j + 1) % vn];
32 if (v0 > v1)
33 std::swap(v0, v1);
34 temp.push_back(std::make_tuple(v0, v1, i, j));
35 }
36 hmi.faces[i].es.resize(vn);
37 }
38 std::sort(temp.begin(), temp.end());
39 hmi.edges.reserve(temp.size() / 2);
40 uint32_t E_num = 0;
41 Edge e;
42 e.boundary = true;
43 e.vs.resize(2);
44 for (uint32_t i = 0; i < temp.size(); ++i)
45 {
46 if (i == 0 || (i != 0 && (std::get<0>(temp[i]) != std::get<0>(temp[i - 1]) || std::get<1>(temp[i]) != std::get<1>(temp[i - 1]))))
47 {
48 e.id = E_num;
49 E_num++;
50 e.vs[0] = std::get<0>(temp[i]);
51 e.vs[1] = std::get<1>(temp[i]);
52 hmi.edges.push_back(e);
53 }
54 else if (i != 0 && (std::get<0>(temp[i]) == std::get<0>(temp[i - 1]) && std::get<1>(temp[i]) == std::get<1>(temp[i - 1])))
55 hmi.edges[E_num - 1].boundary = false;
56
57 hmi.faces[std::get<2>(temp[i])].es[std::get<3>(temp[i])] = E_num - 1;
58 }
59 // boundary
60 for (auto &v : hmi.vertices)
61 v.boundary = false;
62 for (uint32_t i = 0; i < hmi.edges.size(); ++i)
63 if (hmi.edges[i].boundary)
64 {
65 hmi.vertices[hmi.edges[i].vs[0]].boundary = hmi.vertices[hmi.edges[i].vs[1]].boundary = true;
66 }
67 }
68 else if (hmi.type == MeshType::HEX)
69 {
70
71 std::vector<std::vector<uint32_t>> total_fs(hmi.elements.size() * 6);
72 std::vector<std::tuple<uint32_t, uint32_t, uint32_t, uint32_t, uint32_t, uint32_t, uint32_t>> tempF(hmi.elements.size() * 6);
73 std::vector<uint32_t> vs(4);
74 for (uint32_t i = 0; i < hmi.elements.size(); ++i)
75 {
76 for (short j = 0; j < 6; j++)
77 {
78 for (short k = 0; k < 4; k++)
79 vs[k] = hmi.elements[i].vs[hex_face_table[j][k]];
80 uint32_t id = 6 * i + j;
81 total_fs[id] = vs;
82 std::sort(vs.begin(), vs.end());
83 tempF[id] = std::make_tuple(vs[0], vs[1], vs[2], vs[3], id, i, j);
84 }
85 hmi.elements[i].fs.resize(6);
86 }
87 std::sort(tempF.begin(), tempF.end());
88 hmi.faces.reserve(tempF.size() / 3);
89 Face f;
90 f.boundary = true;
91 uint32_t F_num = 0;
92 for (uint32_t i = 0; i < tempF.size(); ++i)
93 {
94 if (i == 0 || (i != 0 && (std::get<0>(tempF[i]) != std::get<0>(tempF[i - 1]) || std::get<1>(tempF[i]) != std::get<1>(tempF[i - 1]) || std::get<2>(tempF[i]) != std::get<2>(tempF[i - 1]) || std::get<3>(tempF[i]) != std::get<3>(tempF[i - 1]))))
95 {
96 f.id = F_num;
97 F_num++;
98 f.vs = total_fs[std::get<4>(tempF[i])];
99 hmi.faces.push_back(f);
100 }
101 else if (i != 0 && (std::get<0>(tempF[i]) == std::get<0>(tempF[i - 1]) && std::get<1>(tempF[i]) == std::get<1>(tempF[i - 1]) && std::get<2>(tempF[i]) == std::get<2>(tempF[i - 1]) && std::get<3>(tempF[i]) == std::get<3>(tempF[i - 1])))
102 hmi.faces[F_num - 1].boundary = false;
103
104 hmi.elements[std::get<5>(tempF[i])].fs[std::get<6>(tempF[i])] = F_num - 1;
105 }
106
107 std::vector<std::tuple<uint32_t, uint32_t, uint32_t, uint32_t>> temp(hmi.faces.size() * 4);
108 for (uint32_t i = 0; i < hmi.faces.size(); ++i)
109 {
110 for (uint32_t j = 0; j < 4; ++j)
111 {
112 uint32_t v0 = hmi.faces[i].vs[j], v1 = hmi.faces[i].vs[(j + 1) % 4];
113 if (v0 > v1)
114 std::swap(v0, v1);
115 temp[4 * i + j] = std::make_tuple(v0, v1, i, j);
116 }
117 hmi.faces[i].es.resize(4);
118 }
119 std::sort(temp.begin(), temp.end());
120 hmi.edges.reserve(temp.size() / 2);
121 uint32_t E_num = 0;
122 Edge e;
123 e.boundary = false;
124 e.vs.resize(2);
125 for (uint32_t i = 0; i < temp.size(); ++i)
126 {
127 if (i == 0 || (i != 0 && (std::get<0>(temp[i]) != std::get<0>(temp[i - 1]) || std::get<1>(temp[i]) != std::get<1>(temp[i - 1]))))
128 {
129 e.id = E_num;
130 E_num++;
131 e.vs[0] = std::get<0>(temp[i]);
132 e.vs[1] = std::get<1>(temp[i]);
133 hmi.edges.push_back(e);
134 }
135 hmi.faces[std::get<2>(temp[i])].es[std::get<3>(temp[i])] = E_num - 1;
136 }
137 // boundary
138 for (auto &v : hmi.vertices)
139 v.boundary = false;
140 for (uint32_t i = 0; i < hmi.faces.size(); ++i)
141 if (hmi.faces[i].boundary)
142 for (uint32_t j = 0; j < 4; ++j)
143 {
144 uint32_t eid = hmi.faces[i].es[j];
145 hmi.edges[eid].boundary = true;
146 hmi.vertices[hmi.edges[eid].vs[0]].boundary = hmi.vertices[hmi.edges[eid].vs[1]].boundary = true;
147 }
148 }
149 else if (hmi.type == MeshType::HYB || hmi.type == MeshType::TET)
150 {
151 vector<bool> bf_flag(hmi.faces.size(), false);
152 for (auto h : hmi.elements)
153 for (auto f : h.fs)
154 bf_flag[f] = !bf_flag[f];
155 for (auto &f : hmi.faces)
156 f.boundary = bf_flag[f.id];
157
158 std::vector<std::tuple<uint32_t, uint32_t, uint32_t, uint32_t>> temp;
159 for (uint32_t i = 0; i < hmi.faces.size(); ++i)
160 {
161 int fl = hmi.faces[i].vs.size();
162 for (uint32_t j = 0; j < hmi.faces[i].vs.size(); ++j)
163 {
164 uint32_t v0 = hmi.faces[i].vs[j], v1 = hmi.faces[i].vs[(j + 1) % fl];
165 if (v0 > v1)
166 std::swap(v0, v1);
167 temp.push_back(std::make_tuple(v0, v1, i, j));
168 }
169 hmi.faces[i].es.resize(fl);
170 }
171 std::sort(temp.begin(), temp.end());
172 hmi.edges.reserve(temp.size() / 2);
173 uint32_t E_num = 0;
174 Edge e;
175 e.boundary = false;
176 e.vs.resize(2);
177 for (uint32_t i = 0; i < temp.size(); ++i)
178 {
179 if (i == 0 || (i != 0 && (std::get<0>(temp[i]) != std::get<0>(temp[i - 1]) || std::get<1>(temp[i]) != std::get<1>(temp[i - 1]))))
180 {
181 e.id = E_num;
182 E_num++;
183 e.vs[0] = std::get<0>(temp[i]);
184 e.vs[1] = std::get<1>(temp[i]);
185 hmi.edges.push_back(e);
186 }
187 hmi.faces[std::get<2>(temp[i])].es[std::get<3>(temp[i])] = E_num - 1;
188 }
189 // boundary
190 for (auto &v : hmi.vertices)
191 v.boundary = false;
192 for (uint32_t i = 0; i < hmi.faces.size(); ++i)
193 if (hmi.faces[i].boundary)
194 for (uint32_t j = 0; j < hmi.faces[i].vs.size(); ++j)
195 {
196 uint32_t eid = hmi.faces[i].es[j];
197 hmi.edges[eid].boundary = true;
198 hmi.vertices[hmi.faces[i].vs[j]].boundary = true;
199 }
200 }
201 // f_nhs;
202 for (auto &f : hmi.faces)
203 f.neighbor_hs.clear();
204 for (uint32_t i = 0; i < hmi.elements.size(); i++)
205 {
206 for (uint32_t j = 0; j < hmi.elements[i].fs.size(); j++)
207 {
208 hmi.faces[hmi.elements[i].fs[j]].neighbor_hs.push_back(i);
209 }
210 }
211 // e_nfs, v_nfs
212 for (auto &e : hmi.edges)
213 e.neighbor_fs.clear();
214 for (auto &v : hmi.vertices)
215 v.neighbor_fs.clear();
216 for (uint32_t i = 0; i < hmi.faces.size(); i++)
217 {
218 for (uint32_t j = 0; j < hmi.faces[i].es.size(); j++)
219 hmi.edges[hmi.faces[i].es[j]].neighbor_fs.push_back(i);
220 for (uint32_t j = 0; j < hmi.faces[i].vs.size(); j++)
221 hmi.vertices[hmi.faces[i].vs[j]].neighbor_fs.push_back(i);
222 }
223 // v_nes, v_nvs
224 for (auto &v : hmi.vertices)
225 {
226 v.neighbor_es.clear();
227 v.neighbor_vs.clear();
228 }
229 for (uint32_t i = 0; i < hmi.edges.size(); i++)
230 {
231 uint32_t v0 = hmi.edges[i].vs[0], v1 = hmi.edges[i].vs[1];
232 hmi.vertices[v0].neighbor_es.push_back(i);
233 hmi.vertices[v1].neighbor_es.push_back(i);
234 hmi.vertices[v0].neighbor_vs.push_back(v1);
235 hmi.vertices[v1].neighbor_vs.push_back(v0);
236 }
237 // e_nhs
238 for (auto &e : hmi.edges)
239 e.neighbor_hs.clear();
240 for (auto &ele : hmi.elements)
241 ele.es.clear();
242 for (uint32_t i = 0; i < hmi.edges.size(); i++)
243 {
244 std::vector<uint32_t> nhs;
245 for (uint32_t j = 0; j < hmi.edges[i].neighbor_fs.size(); j++)
246 {
247 uint32_t nfid = hmi.edges[i].neighbor_fs[j];
248 nhs.insert(nhs.end(), hmi.faces[nfid].neighbor_hs.begin(), hmi.faces[nfid].neighbor_hs.end());
249 }
250 std::sort(nhs.begin(), nhs.end());
251 nhs.erase(std::unique(nhs.begin(), nhs.end()), nhs.end());
252 hmi.edges[i].neighbor_hs = nhs;
253 for (auto nhid : nhs)
254 hmi.elements[nhid].es.push_back(i);
255 }
256 // v_nhs; ordering fs for hex
257 if (hmi.type != MeshType::HYB && hmi.type != MeshType::TET)
258 return;
259
260 for (auto &v : hmi.vertices)
261 v.neighbor_hs.clear();
262 for (uint32_t i = 0; i < hmi.elements.size(); i++)
263 {
264 vector<uint32_t> vs;
265 for (auto fid : hmi.elements[i].fs)
266 vs.insert(vs.end(), hmi.faces[fid].vs.begin(), hmi.faces[fid].vs.end());
267 sort(vs.begin(), vs.end());
268 vs.erase(unique(vs.begin(), vs.end()), vs.end());
269
270 bool degree3 = true;
271 for (auto vid : vs)
272 {
273 int nv = 0;
274 for (auto nvid : hmi.vertices[vid].neighbor_vs)
275 if (find(vs.begin(), vs.end(), nvid) != vs.end())
276 nv++;
277 if (nv != 3)
278 {
279 degree3 = false;
280 break;
281 }
282 }
283
284 if (hmi.elements[i].hex && (vs.size() != 8 || !degree3))
285 hmi.elements[i].hex = false;
286
287 hmi.elements[i].vs.clear();
288
289 if (hmi.elements[i].hex)
290 {
291 int top_fid = hmi.elements[i].fs[0];
292 hmi.elements[i].vs = hmi.faces[top_fid].vs;
293
294 std::set<uint32_t> s_model(vs.begin(), vs.end());
295 std::set<uint32_t> s_pattern(hmi.faces[top_fid].vs.begin(), hmi.faces[top_fid].vs.end());
296 vector<uint32_t> vs_left;
297 std::set_difference(s_model.begin(), s_model.end(), s_pattern.begin(), s_pattern.end(), std::back_inserter(vs_left));
298
299 for (auto vid : hmi.faces[top_fid].vs)
300 for (auto nvid : hmi.vertices[vid].neighbor_vs)
301 if (find(vs_left.begin(), vs_left.end(), nvid) != vs_left.end())
302 {
303 hmi.elements[i].vs.push_back(nvid);
304 break;
305 }
306
307 function<int(vector<uint32_t> &, int &)> WHICH_F = [&](vector<uint32_t> &vs0, int &f_flag) -> int {
308 int which_f = -1;
309 sort(vs0.begin(), vs0.end());
310 bool found_f = false;
311 for (uint32_t j = 0; j < hmi.elements[i].fs.size(); j++)
312 {
313 auto fid = hmi.elements[i].fs[j];
314 vector<uint32_t> vs1 = hmi.faces[fid].vs;
315 sort(vs1.begin(), vs1.end());
316 if (vs0.size() == vs1.size() && std::equal(vs0.begin(), vs0.end(), vs1.begin()))
317 {
318 f_flag = hmi.elements[i].fs_flag[j];
319 which_f = fid;
320 break;
321 }
322 }
323 return which_f;
324 };
325
326 vector<uint32_t> fs;
327 vector<bool> fs_flag;
328 fs_flag.push_back(hmi.elements[i].fs_flag[0]);
329 fs.push_back(top_fid);
330 vector<uint32_t> vs_temp;
331
332 vs_temp.insert(vs_temp.end(), hmi.elements[i].vs.begin() + 4, hmi.elements[i].vs.end());
333 int f_flag = -1;
334 int bottom_fid = WHICH_F(vs_temp, f_flag);
335 fs_flag.push_back(f_flag);
336 fs.push_back(bottom_fid);
337
338 vs_temp.clear();
339 vs_temp.push_back(hmi.elements[i].vs[0]);
340 vs_temp.push_back(hmi.elements[i].vs[1]);
341 vs_temp.push_back(hmi.elements[i].vs[4]);
342 vs_temp.push_back(hmi.elements[i].vs[5]);
343 f_flag = -1;
344 int front_fid = WHICH_F(vs_temp, f_flag);
345 fs_flag.push_back(f_flag);
346 fs.push_back(front_fid);
347
348 vs_temp.clear();
349 vs_temp.push_back(hmi.elements[i].vs[2]);
350 vs_temp.push_back(hmi.elements[i].vs[3]);
351 vs_temp.push_back(hmi.elements[i].vs[6]);
352 vs_temp.push_back(hmi.elements[i].vs[7]);
353 f_flag = -1;
354 int back_fid = WHICH_F(vs_temp, f_flag);
355 fs_flag.push_back(f_flag);
356 fs.push_back(back_fid);
357
358 vs_temp.clear();
359 vs_temp.push_back(hmi.elements[i].vs[1]);
360 vs_temp.push_back(hmi.elements[i].vs[2]);
361 vs_temp.push_back(hmi.elements[i].vs[5]);
362 vs_temp.push_back(hmi.elements[i].vs[6]);
363 f_flag = -1;
364 int left_fid = WHICH_F(vs_temp, f_flag);
365 fs_flag.push_back(f_flag);
366 fs.push_back(left_fid);
367
368 vs_temp.clear();
369 vs_temp.push_back(hmi.elements[i].vs[3]);
370 vs_temp.push_back(hmi.elements[i].vs[0]);
371 vs_temp.push_back(hmi.elements[i].vs[7]);
372 vs_temp.push_back(hmi.elements[i].vs[4]);
373 f_flag = -1;
374 int right_fid = WHICH_F(vs_temp, f_flag);
375 fs_flag.push_back(f_flag);
376 fs.push_back(right_fid);
377
378 hmi.elements[i].fs = fs;
379 hmi.elements[i].fs_flag = fs_flag;
380 }
381 else
382 hmi.elements[i].vs = vs;
383 for (uint32_t j = 0; j < hmi.elements[i].vs.size(); j++)
384 hmi.vertices[hmi.elements[i].vs[j]].neighbor_hs.push_back(i);
385 }
386 // matrix representation of tet mesh
387 if (hmi.type == MeshType::TET)
388 {
389 hmi.EV.resize(2, hmi.edges.size());
390 for (const auto &e : hmi.edges)
391 {
392 hmi.EV(0, e.id) = e.vs[0];
393 hmi.EV(1, e.id) = e.vs[1];
394 }
395 hmi.FV.resize(3, hmi.faces.size());
396 hmi.FE.resize(3, hmi.faces.size());
397 hmi.FH.resize(2, hmi.faces.size());
398 hmi.FHi.resize(2, hmi.faces.size());
399 for (const auto &f : hmi.faces)
400 {
401 hmi.FV(0, f.id) = f.vs[0];
402 hmi.FV(1, f.id) = f.vs[1];
403 hmi.FV(2, f.id) = f.vs[2];
404
405 hmi.FE(0, f.id) = f.es[0];
406 hmi.FE(1, f.id) = f.es[1];
407 hmi.FE(2, f.id) = f.es[2];
408
409 hmi.FH(0, f.id) = f.neighbor_hs[0];
410 for (int i = 0; i < hmi.elements[f.neighbor_hs[0]].fs.size(); i++)
411 if (f.id == hmi.elements[f.neighbor_hs[0]].fs[i])
412 hmi.FHi(0, f.id) = i;
413
414 hmi.FH(1, f.id) = -1;
415 hmi.FHi(1, f.id) = -1;
416 if (f.neighbor_hs.size() == 2)
417 {
418 hmi.FH(1, f.id) = f.neighbor_hs[1];
419 for (int i = 0; i < hmi.elements[f.neighbor_hs[1]].fs.size(); i++)
420 if (f.id == hmi.elements[f.neighbor_hs[1]].fs[i])
421 hmi.FHi(1, f.id) = i;
422 }
423 }
424 hmi.HV.resize(4, hmi.elements.size());
425 hmi.HF.resize(4, hmi.elements.size());
426 for (const auto &h : hmi.elements)
427 {
428 hmi.HV(0, h.id) = h.vs[0];
429 hmi.HV(1, h.id) = h.vs[1];
430 hmi.HV(2, h.id) = h.vs[2];
431 hmi.HV(3, h.id) = h.vs[3];
432
433 hmi.HF(0, h.id) = h.fs[0];
434 hmi.HF(1, h.id) = h.fs[1];
435 hmi.HF(2, h.id) = h.fs[2];
436 hmi.HF(3, h.id) = h.fs[3];
437 }
438 }
439
440 // boundary flags for hybrid mesh
441 std::vector<bool> bv_flag(hmi.vertices.size(), false), be_flag(hmi.edges.size(), false), bf_flag(hmi.faces.size(), false);
442 for (auto f : hmi.faces)
443 if (f.boundary && hmi.elements[f.neighbor_hs[0]].hex)
444 bf_flag[f.id] = true;
445 else if (!f.boundary)
446 {
447 int ele0 = f.neighbor_hs[0], ele1 = f.neighbor_hs[1];
448 if ((hmi.elements[ele0].hex && !hmi.elements[ele1].hex) || (!hmi.elements[ele0].hex && hmi.elements[ele1].hex))
449 bf_flag[f.id] = true;
450 }
451 for (uint32_t i = 0; i < hmi.faces.size(); ++i)
452 if (bf_flag[i])
453 for (uint32_t j = 0; j < hmi.faces[i].vs.size(); ++j)
454 {
455 uint32_t eid = hmi.faces[i].es[j];
456 be_flag[eid] = true;
457 bv_flag[hmi.faces[i].vs[j]] = true;
458 }
459 // boundary_hex for hybrid mesh
460 for (auto &v : hmi.vertices)
461 v.boundary_hex = false;
462 for (auto &e : hmi.edges)
463 e.boundary_hex = false;
464 for (auto &f : hmi.faces)
465 f.boundary_hex = bf_flag[f.id];
466 for (auto &ele : hmi.elements)
467 if (ele.hex)
468 {
469 for (auto vid : ele.vs)
470 {
471 if (!bv_flag[vid])
472 continue;
473 int fn = 0;
474 for (auto nfid : hmi.vertices[vid].neighbor_fs)
475 if (bf_flag[nfid] && std::find(ele.fs.begin(), ele.fs.end(), nfid) != ele.fs.end())
476 fn++;
477 if (fn == 3)
478 hmi.vertices[vid].boundary_hex = true;
479 }
480 for (auto eid : ele.es)
481 {
482 if (!be_flag[eid])
483 continue;
484 int fn = 0;
485 for (auto nfid : hmi.edges[eid].neighbor_fs)
486 if (bf_flag[nfid] && std::find(ele.fs.begin(), ele.fs.end(), nfid) != ele.fs.end())
487 fn++;
488 if (fn == 2)
489 hmi.edges[eid].boundary_hex = true;
490 }
491 }
492}
493void MeshProcessing3D::reorder_hex_mesh_propogation(Mesh3DStorage &hmi)
494{
495 // connected components
496 vector<bool> H_tag(hmi.elements.size(), false);
497 vector<vector<uint32_t>> Groups;
498 while (true)
499 {
500 vector<uint32_t> group;
501 int shid = -1;
502 for (uint32_t i = 0; i < H_tag.size(); i++)
503 if (!H_tag[i])
504 {
505 shid = i;
506 break;
507 }
508 if (shid == -1)
509 break;
510 group.push_back(shid);
511 H_tag[shid] = true;
512 vector<uint32_t> group_ = group;
513 while (true)
514 {
515 vector<uint32_t> pool;
516 for (auto hid : group_)
517 {
518 vector<vector<uint32_t>> Fvs(6), fvs_sorted;
519 for (uint32_t i = 0; i < 6; i++)
520 for (uint32_t j = 0; j < 4; j++)
521 Fvs[i].push_back(hmi.elements[hid].vs[hex_face_table[i][j]]);
522 fvs_sorted = Fvs;
523 for (auto &vs : fvs_sorted)
524 sort(vs.begin(), vs.end());
525
526 for (auto fid : hmi.elements[hid].fs)
527 if (!hmi.faces[fid].boundary)
528 {
529 int nhid = hmi.faces[fid].neighbor_hs[0];
530 if (nhid == hid)
531 nhid = hmi.faces[fid].neighbor_hs[1];
532
533 if (!H_tag[nhid])
534 {
535 pool.push_back(nhid);
536 H_tag[nhid] = true;
537
538 vector<uint32_t> fvs = hmi.faces[fid].vs;
539 sort(fvs.begin(), fvs.end());
540
541 int f_ind = -1;
542 for (uint32_t i = 0; i < 6; i++)
543 if (std::equal(fvs.begin(), fvs.end(), fvs_sorted[i].begin()))
544 {
545 f_ind = i;
546 break;
547 }
548 vector<uint32_t> hvs = hmi.elements[nhid].vs;
549 hmi.elements[nhid].vs.clear();
550 vector<uint32_t> topvs = Fvs[f_ind];
551 std::reverse(topvs.begin(), topvs.end());
552 vector<uint32_t> bottomvs;
553 for (uint32_t i = 0; i < 4; i++)
554 {
555 for (auto nvid : hmi.vertices[topvs[i]].neighbor_vs)
556 if (nvid != topvs[(i + 3) % 4] && nvid != topvs[(i + 1) % 4]
557 && std::find(hvs.begin(), hvs.end(), nvid) != hvs.end())
558 {
559 bottomvs.push_back(nvid);
560 break;
561 }
562 }
563 hmi.elements[nhid].vs = topvs;
564 hmi.elements[nhid].vs.insert(hmi.elements[nhid].vs.end(), bottomvs.begin(), bottomvs.end());
565 }
566 }
567 }
568 if (pool.size())
569 {
570 group_ = pool;
571 group.insert(group.end(), pool.begin(), pool.end());
572 pool.size();
573 }
574 else
575 break;
576 }
577 Groups.push_back(group);
578 }
579 // direction
580 for (auto group : Groups)
581 {
582 Mesh_Quality mq1, mq2;
583 Mesh3DStorage m1, m2;
584 m2.type = m1.type = MeshType::HEX;
585
586 m2.points = m1.points = hmi.points;
587 for (auto hid : group)
588 m1.elements.push_back(hmi.elements[hid]);
589 scaled_jacobian(m1, mq1);
590 if (mq1.min_Jacobian > 0)
591 continue;
592 m2.elements = m1.elements;
593 for (auto &h : m2.elements)
594 {
595 swap(h.vs[1], h.vs[3]);
596 swap(h.vs[5], h.vs[7]);
597 }
598 scaled_jacobian(m2, mq2);
599 if (mq2.ave_Jacobian > mq1.ave_Jacobian)
600 {
601 for (auto &h : m2.elements)
602 hmi.elements[h.id] = h;
603 }
604 }
605}
606bool MeshProcessing3D::scaled_jacobian(Mesh3DStorage &hmi, Mesh_Quality &mq)
607{
608 if (hmi.type != MeshType::HEX)
609 return false;
610
611 mq.ave_Jacobian = 0;
612 mq.min_Jacobian = 1;
613 mq.deviation_Jacobian = 0;
614 mq.V_Js.resize(hmi.elements.size() * 8);
615 mq.V_Js.setZero();
616 mq.H_Js.resize(hmi.elements.size());
617 mq.H_Js.setZero();
618
619 for (uint32_t i = 0; i < hmi.elements.size(); i++)
620 {
621 double hex_minJ = 1;
622 for (uint32_t j = 0; j < 8; j++)
623 {
624 uint32_t v0, v1, v2, v3;
625 v0 = hex_tetra_table[j][0];
626 v1 = hex_tetra_table[j][1];
627 v2 = hex_tetra_table[j][2];
628 v3 = hex_tetra_table[j][3];
629
630 Vector3d c0 = hmi.points.col(hmi.elements[i].vs[v0]);
631 Vector3d c1 = hmi.points.col(hmi.elements[i].vs[v1]);
632 Vector3d c2 = hmi.points.col(hmi.elements[i].vs[v2]);
633 Vector3d c3 = hmi.points.col(hmi.elements[i].vs[v3]);
634
635 double jacobian_value = a_jacobian(c0, c1, c2, c3);
636
637 if (hex_minJ > jacobian_value)
638 hex_minJ = jacobian_value;
639
640 uint32_t id = 8 * i + j;
641 mq.V_Js[id] = jacobian_value;
642 }
643 mq.H_Js[i] = hex_minJ;
644 mq.ave_Jacobian += hex_minJ;
645 if (mq.min_Jacobian > hex_minJ)
646 mq.min_Jacobian = hex_minJ;
647 }
648 mq.ave_Jacobian /= hmi.elements.size();
649 for (int i = 0; i < mq.H_Js.size(); i++)
650 mq.deviation_Jacobian += (mq.H_Js[i] - mq.ave_Jacobian) * (mq.H_Js[i] - mq.ave_Jacobian);
651 mq.deviation_Jacobian /= hmi.elements.size();
652
653 return true;
654}
655double MeshProcessing3D::a_jacobian(Vector3d &v0, Vector3d &v1, Vector3d &v2, Vector3d &v3)
656{
657 Matrix3d Jacobian;
658
659 Jacobian.col(0) = v1 - v0;
660 Jacobian.col(1) = v2 - v0;
661 Jacobian.col(2) = v3 - v0;
662
663 double norm1 = Jacobian.col(0).norm();
664 double norm2 = Jacobian.col(1).norm();
665 double norm3 = Jacobian.col(2).norm();
666
667 double scaled_jacobian = Jacobian.determinant();
668 if (std::abs(norm1) < Jacobian_Precision || std::abs(norm2) < Jacobian_Precision || std::abs(norm3) < Jacobian_Precision)
669 {
670 return scaled_jacobian;
671 }
672 scaled_jacobian /= norm1 * norm2 * norm3;
673 return scaled_jacobian;
674}
675
676void MeshProcessing3D::global_orientation_hexes(Mesh3DStorage &hmi)
677{
678 Mesh3DStorage mesh;
679 mesh.type = MeshType::HEX;
680 mesh.points = hmi.points;
681 mesh.vertices.resize(hmi.vertices.size());
682 for (const auto &v : hmi.vertices)
683 {
684 Vertex v_;
685 v_.id = v.id;
686 v_.v = v.v;
687 mesh.vertices[v.id] = v;
688 }
689 vector<int> Ele_map(hmi.elements.size(), -1), Ele_map_reverse;
690 for (const auto &ele : hmi.elements)
691 {
692 if (!ele.hex)
693 continue;
694 Element ele_;
695 ele_.id = mesh.elements.size();
696 ele_.vs = ele.vs;
697 for (auto vid : ele_.vs)
698 mesh.vertices[vid].neighbor_hs.push_back(ele_.id);
699 mesh.elements.push_back(ele_);
700
701 Ele_map[ele.id] = ele_.id;
702 Ele_map_reverse.push_back(ele.id);
703 }
704
705 build_connectivity(mesh);
706 reorder_hex_mesh_propogation(mesh);
707
708 for (const auto &ele : mesh.elements)
709 hmi.elements[Ele_map_reverse[ele.id]].vs = ele.vs;
710}
711void MeshProcessing3D::refine_catmul_clark_polar(Mesh3DStorage &M, int iter, bool reverse, std::vector<int> &Parents)
712{
713
714 for (int i = 0; i < iter; i++)
715 {
716
717 std::vector<int> Refinement_Levels;
718 ele_subdivison_levels(M, Refinement_Levels);
719
720 Mesh3DStorage M_;
721 M_.type = MeshType::HYB;
722
723 vector<int> E2V(M.edges.size()), F2V(M.faces.size()), Ele2V(M.elements.size());
724
725 int vn = 0;
726 for (auto v : M.vertices)
727 {
728 Vertex v_;
729 v_.id = vn++;
730 v_.v.resize(3);
731 for (int j = 0; j < 3; j++)
732 v_.v[j] = M.points(j, v.id);
733 M_.vertices.push_back(v_);
734 }
735
736 for (auto e : M.edges)
737 {
738 Vertex v;
739 v.id = vn++;
740 v.v.resize(3);
741
742 Vector3d center;
743 center.setZero();
744 for (auto vid : e.vs)
745 center += M.points.col(vid);
746 center /= e.vs.size();
747 for (int j = 0; j < 3; j++)
748 v.v[j] = center[j];
749
750 M_.vertices.push_back(v);
751 E2V[e.id] = v.id;
752 }
753 for (auto f : M.faces)
754 {
755 Vertex v;
756 v.id = vn++;
757 v.v.resize(3);
758
759 Vector3d center;
760 center.setZero();
761 for (auto vid : f.vs)
762 center += M.points.col(vid);
763 center /= f.vs.size();
764 for (int j = 0; j < 3; j++)
765 v.v[j] = center[j];
766
767 M_.vertices.push_back(v);
768 F2V[f.id] = v.id;
769 }
770 for (auto ele : M.elements)
771 {
772 if (!ele.hex)
773 continue;
774 Vertex v;
775 v.id = vn++;
776 v.v = ele.v_in_Kernel;
777
778 M_.vertices.push_back(v);
779 Ele2V[ele.id] = v.id;
780 }
781 // new elements
782 std::vector<std::vector<uint32_t>> total_fs;
783 total_fs.reserve(M.elements.size() * 8 * 6);
784 std::vector<std::tuple<uint32_t, uint32_t, uint32_t, uint32_t, uint32_t, uint32_t, uint32_t>> tempF;
785 tempF.reserve(M.elements.size() * 8 * 6);
786 std::vector<uint32_t> vs(4);
787
788 int elen = 0, fn = 0;
789 for (auto ele : M.elements)
790 {
791 if (ele.hex)
792 {
793 for (auto vid : ele.vs)
794 {
795 // top 4 vs
796 vector<int> top_vs(4);
797 top_vs[0] = vid;
798 int fid = -1;
799 for (auto nfid : M.vertices[vid].neighbor_fs)
800 if (find(ele.fs.begin(), ele.fs.end(), nfid) != ele.fs.end())
801 {
802 fid = nfid;
803 break;
804 }
805 assert(fid != -1);
806 top_vs[2] = F2V[fid];
807
808 int v_ind = find(M.faces[fid].vs.begin(), M.faces[fid].vs.end(), vid) - M.faces[fid].vs.begin();
809 int e_pre = M.faces[fid].es[(v_ind - 1 + 4) % 4];
810 int e_aft = M.faces[fid].es[v_ind];
811 top_vs[1] = E2V[e_pre];
812 top_vs[3] = E2V[e_aft];
813 // bottom 4 vs
814 vector<int> bottom_vs(4);
815
816 auto nvs = M.vertices[vid].neighbor_vs;
817 int e_per = -1;
818 for (auto nvid : nvs)
819 if (find(ele.vs.begin(), ele.vs.end(), nvid) != ele.vs.end())
820 {
821 if (nvid != M.edges[e_pre].vs[0] && nvid != M.edges[e_pre].vs[1] && nvid != M.edges[e_aft].vs[0] && nvid != M.edges[e_aft].vs[1])
822 {
823 vector<uint32_t> sharedes, es0 = M.vertices[vid].neighbor_es, es1 = M.vertices[nvid].neighbor_es;
824 sort(es0.begin(), es0.end());
825 sort(es1.begin(), es1.end());
826 set_intersection(es0.begin(), es0.end(), es1.begin(), es1.end(), back_inserter(sharedes));
827 assert(sharedes.size());
828 e_per = sharedes[0];
829 break;
830 }
831 }
832
833 assert(e_per != -1);
834 bottom_vs[0] = E2V[e_per];
835 bottom_vs[2] = Ele2V[ele.id];
836
837 int f_pre = -1;
838 vector<uint32_t> sharedfs, fs0 = M.edges[e_pre].neighbor_fs, fs1 = M.edges[e_per].neighbor_fs;
839 sort(fs0.begin(), fs0.end());
840 sort(fs1.begin(), fs1.end());
841 set_intersection(fs0.begin(), fs0.end(), fs1.begin(), fs1.end(), back_inserter(sharedfs));
842 for (auto sfid : sharedfs)
843 if (find(ele.fs.begin(), ele.fs.end(), sfid) != ele.fs.end())
844 {
845 f_pre = sfid;
846 break;
847 }
848 assert(f_pre != -1);
849
850 int f_aft = -1;
851 sharedfs.clear();
852 fs0 = M.edges[e_aft].neighbor_fs;
853 fs1 = M.edges[e_per].neighbor_fs;
854 sort(fs0.begin(), fs0.end());
855 sort(fs1.begin(), fs1.end());
856 set_intersection(fs0.begin(), fs0.end(), fs1.begin(), fs1.end(), back_inserter(sharedfs));
857 for (auto sfid : sharedfs)
858 if (find(ele.fs.begin(), ele.fs.end(), sfid) != ele.fs.end())
859 {
860 f_aft = sfid;
861 break;
862 }
863 assert(f_aft != -1);
864
865 bottom_vs[1] = F2V[f_pre];
866 bottom_vs[3] = F2V[f_aft];
867
868 vector<int> ele_vs = top_vs;
869 ele_vs.insert(ele_vs.end(), bottom_vs.begin(), bottom_vs.end());
870 // fs
871 for (short j = 0; j < 6; j++)
872 {
873 for (short k = 0; k < 4; k++)
874 vs[k] = ele_vs[hex_face_table[j][k]];
875 total_fs.push_back(vs);
876 std::sort(vs.begin(), vs.end());
877 tempF.push_back(std::make_tuple(vs[0], vs[1], vs[2], vs[3], fn++, elen, j));
878 }
879 // new ele
880 Element ele_;
881 ele_.id = elen++;
882 ele_.fs.resize(6, -1);
883 ele_.fs_flag.resize(6, 1);
884
885 ele_.hex = true;
886 ele_.v_in_Kernel.resize(3);
887
888 Vector3d center;
889 center.setZero();
890 for (auto evid : ele_vs)
891 for (int j = 0; j < 3; j++)
892 center[j] += M_.vertices[evid].v[j];
893 center /= ele_vs.size();
894 for (int j = 0; j < 3; j++)
895 ele_.v_in_Kernel[j] = center[j];
896
897 M_.elements.push_back(ele_);
898 Parents.push_back(ele.id);
899 }
900 }
901 else
902 {
903 int level = Refinement_Levels[ele.id];
904 if (reverse)
905 level = 1;
906 // local_V2V
907 std::vector<std::vector<int>> local_V2Vs;
908 std::map<int, int> local_vi_map;
909 for (auto vid : ele.vs)
910 {
911 std::vector<int> v2v;
912 v2v.push_back(vid);
913 for (int r = 0; r < level; r++)
914 {
915 Vertex v_;
916 v_.id = vn++;
917 v_.v.resize(3);
918
919 for (int j = 0; j < 3; j++)
920 v_.v[j] = M_.vertices[vid].v[j] + (ele.v_in_Kernel[j] - M_.vertices[vid].v[j]) * (r + 1.0) / (double)(level + 1);
921 if (reverse)
922 {
923 for (int j = 0; j < 3; j++)
924 v_.v[j] = M_.vertices[vid].v[j] + (M_.vertices[vid].v[j] - ele.v_in_Kernel[j]) * (r + 1.0) / (double)(level + 1);
925 }
926 M_.vertices.push_back(v_);
927 v2v.push_back(v_.id);
928 }
929 local_vi_map[vid] = local_V2Vs.size();
930 local_V2Vs.push_back(v2v);
931 }
932 // local_E2V
933 vector<uint32_t> es;
934 for (auto fid : ele.fs)
935 es.insert(es.end(), M.faces[fid].es.begin(), M.faces[fid].es.end());
936 sort(es.begin(), es.end());
937 es.erase(unique(es.begin(), es.end()), es.end());
938
939 std::vector<std::vector<int>> local_E2Vs;
940 std::map<int, int> local_ei_map;
941 for (auto eid : es)
942 {
943 std::vector<int> e2v;
944 e2v.push_back(E2V[eid]);
945 for (int r = 0; r < level; r++)
946 {
947 Vertex v_;
948 v_.id = vn++;
949 v_.v.resize(3);
950
951 Vector3d center;
952 center.setZero();
953 for (auto vid : M.edges[eid].vs)
954 for (int j = 0; j < 3; j++)
955 center[j] += M_.vertices[local_V2Vs[local_vi_map[vid]][r + 1]].v[j];
956 center /= M.edges[eid].vs.size();
957 for (int j = 0; j < 3; j++)
958 v_.v[j] = center[j];
959
960 M_.vertices.push_back(v_);
961 e2v.push_back(v_.id);
962 }
963 local_ei_map[eid] = local_E2Vs.size();
964 local_E2Vs.push_back(e2v);
965 }
966 // local_F2V
967 std::vector<std::vector<int>> local_F2Vs;
968 std::map<int, int> local_fi_map;
969 for (auto fid : ele.fs)
970 {
971 std::vector<int> f2v;
972 f2v.push_back(F2V[fid]);
973 for (int r = 0; r < level; r++)
974 {
975 Vertex v_;
976 v_.id = vn++;
977 v_.v.resize(3);
978
979 Vector3d center;
980 center.setZero();
981 for (auto vid : M.faces[fid].vs)
982 for (int j = 0; j < 3; j++)
983 center[j] += M_.vertices[local_V2Vs[local_vi_map[vid]][r + 1]].v[j];
984 center /= M.faces[fid].vs.size();
985 for (int j = 0; j < 3; j++)
986 v_.v[j] = center[j];
987
988 M_.vertices.push_back(v_);
989 f2v.push_back(v_.id);
990 }
991 local_fi_map[fid] = local_F2Vs.size();
992 local_F2Vs.push_back(f2v);
993 }
994 // polyhedron fs
995 int local_fn = 0;
996 for (auto fid : ele.fs)
997 {
998 auto &fvs = M.faces[fid].vs;
999 auto &fes = M.faces[fid].es;
1000 int fvn = M.faces[fid].vs.size();
1001 for (uint32_t j = 0; j < fvn; j++)
1002 {
1003 vs[0] = local_E2Vs[local_ei_map[fes[(j - 1 + fvn) % fvn]]][level];
1004 vs[1] = local_V2Vs[local_vi_map[fvs[j]]][level];
1005 vs[2] = local_E2Vs[local_ei_map[fes[j]]][level];
1006 vs[3] = local_F2Vs[local_fi_map[fid]][level];
1007
1008 if (reverse)
1009 {
1010 vs[0] = local_E2Vs[local_ei_map[fes[(j - 1 + fvn) % fvn]]][0];
1011 vs[1] = local_V2Vs[local_vi_map[fvs[j]]][0];
1012 vs[2] = local_E2Vs[local_ei_map[fes[j]]][0];
1013 vs[3] = local_F2Vs[local_fi_map[fid]][0];
1014 }
1015
1016 total_fs.push_back(vs);
1017 std::sort(vs.begin(), vs.end());
1018 tempF.push_back(std::make_tuple(vs[0], vs[1], vs[2], vs[3], fn++, elen, local_fn++));
1019 }
1020 }
1021 // polyhedron
1022 Element ele_;
1023 ele_.id = elen++;
1024 ele_.fs.resize(local_fn, -1);
1025 ele_.fs_flag.resize(local_fn, 1);
1026
1027 ele_.hex = false;
1028 ele_.v_in_Kernel = ele.v_in_Kernel;
1029 M_.elements.push_back(ele_);
1030 Parents.push_back(ele.id);
1031 // hex
1032 for (int r = 0; r < level; r++)
1033 {
1034 for (auto fid : ele.fs)
1035 {
1036 auto &fvs = M.faces[fid].vs;
1037 auto &fes = M.faces[fid].es;
1038 int fvn = M.faces[fid].vs.size();
1039 for (uint32_t j = 0; j < fvn; j++)
1040 {
1041 vector<int> ele_vs(8);
1042 ele_vs[0] = local_E2Vs[local_ei_map[fes[(j - 1 + fvn) % fvn]]][r + 1];
1043 ele_vs[1] = local_V2Vs[local_vi_map[fvs[j]]][r + 1];
1044 ele_vs[2] = local_E2Vs[local_ei_map[fes[j]]][r + 1];
1045 ele_vs[3] = local_F2Vs[local_fi_map[fid]][r + 1];
1046
1047 ele_vs[4] = local_E2Vs[local_ei_map[fes[(j - 1 + fvn) % fvn]]][r];
1048 ele_vs[5] = local_V2Vs[local_vi_map[fvs[j]]][r];
1049 ele_vs[6] = local_E2Vs[local_ei_map[fes[j]]][r];
1050 ele_vs[7] = local_F2Vs[local_fi_map[fid]][r];
1051
1052 // fs
1053 for (short j = 0; j < 6; j++)
1054 {
1055 for (short k = 0; k < 4; k++)
1056 vs[k] = ele_vs[hex_face_table[j][k]];
1057 total_fs.push_back(vs);
1058 std::sort(vs.begin(), vs.end());
1059 tempF.push_back(std::make_tuple(vs[0], vs[1], vs[2], vs[3], fn++, elen, j));
1060 }
1061 // hex
1062 Element ele_;
1063 ele_.id = elen++;
1064 ele_.fs.resize(6, -1);
1065 ele_.fs_flag.resize(6, 1);
1066 ele_.hex = true;
1067 ele_.v_in_Kernel.resize(3);
1068
1069 Vector3d center;
1070 center.setZero();
1071 for (auto vid : ele_vs)
1072 for (int j = 0; j < 3; j++)
1073 center[j] += M_.vertices[vid].v[j];
1074 center /= ele_vs.size();
1075 for (int j = 0; j < 3; j++)
1076 ele_.v_in_Kernel[j] = center[j];
1077 M_.elements.push_back(ele_);
1078 Parents.push_back(ele.id);
1079 }
1080 }
1081 }
1082 }
1083 }
1084 // Fs
1085 std::sort(tempF.begin(), tempF.end());
1086 M_.faces.reserve(tempF.size() / 3);
1087 Face f;
1088 f.boundary = true;
1089 uint32_t F_num = 0;
1090 for (uint32_t i = 0; i < tempF.size(); ++i)
1091 {
1092 if (i == 0 || (i != 0 && (std::get<0>(tempF[i]) != std::get<0>(tempF[i - 1]) || std::get<1>(tempF[i]) != std::get<1>(tempF[i - 1]) || std::get<2>(tempF[i]) != std::get<2>(tempF[i - 1]) || std::get<3>(tempF[i]) != std::get<3>(tempF[i - 1]))))
1093 {
1094 f.id = F_num;
1095 F_num++;
1096 f.vs = total_fs[std::get<4>(tempF[i])];
1097 M_.faces.push_back(f);
1098 }
1099 else if (i != 0 && (std::get<0>(tempF[i]) == std::get<0>(tempF[i - 1]) && std::get<1>(tempF[i]) == std::get<1>(tempF[i - 1]) && std::get<2>(tempF[i]) == std::get<2>(tempF[i - 1]) && std::get<3>(tempF[i]) == std::get<3>(tempF[i - 1])))
1100 M_.faces[F_num - 1].boundary = false;
1101
1102 M_.elements[std::get<5>(tempF[i])].fs[std::get<6>(tempF[i])] = F_num - 1;
1103 }
1104
1105 M_.points.resize(3, M_.vertices.size());
1106 for (auto v : M_.vertices)
1107 for (int j = 0; j < 3; j++)
1108 M_.points(j, v.id) = v.v[j];
1109
1110 build_connectivity(M_);
1111 orient_volume_mesh(M_);
1112 build_connectivity(M_);
1113
1114 M = M_;
1115 }
1116}
1117void MeshProcessing3D::refine_red_refinement_tet(Mesh3DStorage &M, int iter)
1118{
1119 for (int i = 0; i < iter; i++)
1120 {
1121
1122 Mesh3DStorage M_;
1123 M_.type = MeshType::TET;
1124
1125 vector<int> E2V(M.edges.size());
1126
1127 int vn = 0;
1128 for (const auto &v : M.vertices)
1129 {
1130 Vertex v_;
1131 v_.id = vn++;
1132 v_.v.resize(3);
1133 for (int j = 0; j < 3; j++)
1134 v_.v[j] = M.points(j, v.id);
1135 M_.vertices.push_back(v_);
1136 }
1137
1138 for (const auto &e : M.edges)
1139 {
1140 Vertex v;
1141 v.id = vn++;
1142 v.v.resize(3);
1143
1144 Vector3d center;
1145 center.setZero();
1146 for (auto vid : e.vs)
1147 center += M.points.col(vid);
1148 center /= e.vs.size();
1149 for (int j = 0; j < 3; j++)
1150 v.v[j] = center[j];
1151
1152 M_.vertices.push_back(v);
1153 E2V[e.id] = v.id;
1154 }
1155 // for (const auto &ele : M.elements) {
1156 // Element ele_;
1157 // ele_.id = M_.elements.size();
1158 // ele_.fs.resize(4, -1);
1159 // ele_.fs_flag.resize(4, 1);
1160 // ele_.vs = ele.vs;
1161
1162 // ele_.hex = false;
1163 // ele_.v_in_Kernel.resize(3);
1164
1165 // Vector3d center;
1166 // center.setZero();
1167 // for (const auto &evid : ele_.vs) for (int j = 0; j < 3; j++)center[j] += M_.vertices[evid].v[j];
1168 // center /= ele_.vs.size();
1169 // for (int j = 0; j < 3; j++)ele_.v_in_Kernel[j] = center[j];
1170
1171 // M_.elements.push_back(ele_);
1172 // }
1173
1174 auto shared_edge = [&](int v0, int v1, int &e) -> bool {
1175 auto &es0 = M.vertices[v0].neighbor_es, &es1 = M.vertices[v1].neighbor_es;
1176 std::sort(es0.begin(), es0.end());
1177 std::sort(es1.begin(), es1.end());
1178 std::vector<uint32_t> es;
1179 std::set_intersection(es0.begin(), es0.end(), es1.begin(), es1.end(), back_inserter(es));
1180 if (es.size())
1181 {
1182 e = es[0];
1183 return true;
1184 }
1185 return false;
1186 };
1187
1188 std::vector<bool> e_flag(M.edges.size(), false);
1189
1190 for (const auto &ele : M.elements)
1191 { // 1 --> 8
1192
1193 for (short i = 0; i < 4; i++)
1194 { // four corners
1195 Element ele_;
1196 ele_.id = M_.elements.size();
1197 ele_.fs.resize(4, -1);
1198 ele_.fs_flag.resize(4, 1);
1199
1200 ele_.hex = false;
1201 ele_.v_in_Kernel.resize(3);
1202
1203 ele_.vs.push_back(ele.vs[i]);
1204 for (short j = 0; j < 4; j++)
1205 {
1206 if (j == i)
1207 continue;
1208 int v0 = ele.vs[i], v1 = ele.vs[j];
1209 int e = -1;
1210 if (shared_edge(v0, v1, e))
1211 ele_.vs.push_back(E2V[e]);
1212 }
1213 Vector3d center;
1214 center.setZero();
1215 for (const auto &evid : ele_.vs)
1216 for (int j = 0; j < 3; j++)
1217 center[j] += M_.vertices[evid].v[j];
1218 center /= ele_.vs.size();
1219 for (int j = 0; j < 3; j++)
1220 ele_.v_in_Kernel[j] = center[j];
1221
1222 M_.elements.push_back(ele_);
1223 }
1224
1225 // 6 edges
1226 std::vector<int> edges(6);
1227 for (int i = 0; i < 6; i++)
1228 {
1229 int v0 = ele.vs[tet_edges[i][0]], v1 = ele.vs[tet_edges[i][1]];
1230 int e = -1;
1231 if (shared_edge(v0, v1, e))
1232 edges[i] = e;
1233 }
1234 // cut the inner octahedron along its shortest diagonal, which joins the midpoints
1235 // of two opposite edges (k, 5 - k). A fixed choice of diagonal makes the tets
1236 // degenerate under repeated refinement, since build_connectivity sorts ele.vs.
1237 int diag = 0;
1238 double min_length = std::numeric_limits<double>::max();
1239 for (int k = 0; k < 3; ++k)
1240 {
1241 const double length = (Map<const Vector3d>(M_.vertices[E2V[edges[k]]].v.data()) - Map<const Vector3d>(M_.vertices[E2V[edges[5 - k]]].v.data())).squaredNorm();
1242 if (length < min_length)
1243 {
1244 min_length = length;
1245 diag = k;
1246 }
1247 }
1248 int lv0 = E2V[edges[diag]], lv1 = E2V[edges[5 - diag]];
1249 e_flag[edges[diag]] = true;
1250 e_flag[edges[5 - diag]] = true;
1251 for (short i = 0; i < 4; i++)
1252 { // four faces
1253 Element ele_;
1254 ele_.id = M_.elements.size();
1255 ele_.fs.resize(4, -1);
1256 ele_.fs_flag.resize(4, 1);
1257
1258 ele_.hex = false;
1259 ele_.v_in_Kernel.resize(3);
1260
1261 ele_.vs.push_back(lv0);
1262 ele_.vs.push_back(lv1);
1263 for (short j = 0; j < 3; j++)
1264 {
1265 int c_e = M.faces[ele.fs[i]].es[j];
1266 if (e_flag[c_e])
1267 continue;
1268 ele_.vs.push_back(E2V[c_e]);
1269 }
1270 Vector3d center;
1271 center.setZero();
1272 for (const auto &evid : ele_.vs)
1273 for (int j = 0; j < 3; j++)
1274 center[j] += M_.vertices[evid].v[j];
1275 center /= ele_.vs.size();
1276 for (int j = 0; j < 3; j++)
1277 ele_.v_in_Kernel[j] = center[j];
1278
1279 M_.elements.push_back(ele_);
1280 }
1281 e_flag[edges[diag]] = e_flag[edges[5 - diag]] = false;
1282 }
1283
1284 M_.points.resize(3, M_.vertices.size());
1285 for (auto v : M_.vertices)
1286 for (int j = 0; j < 3; j++)
1287 M_.points(j, v.id) = v.v[j];
1288
1289 // orient tets
1290 const auto &t = M.elements[0].vs;
1291 Vector3d c0 = M.points.col(t[0]);
1292 Vector3d c1 = M.points.col(t[1]);
1293 Vector3d c2 = M.points.col(t[2]);
1294 Vector3d c3 = M.points.col(t[3]);
1295 bool signed_volume = a_jacobian(c0, c1, c2, c3) > 0 ? true : false;
1296
1297 for (auto &ele : M_.elements)
1298 {
1299 c0 = M_.points.col(ele.vs[0]);
1300 c1 = M_.points.col(ele.vs[1]);
1301 c2 = M_.points.col(ele.vs[2]);
1302 c3 = M_.points.col(ele.vs[3]);
1303 bool sign = a_jacobian(c0, c1, c2, c3) > 0 ? true : false;
1304 if (sign != signed_volume)
1305 std::swap(ele.vs[1], ele.vs[3]);
1306 }
1307
1308 // Fs
1309 std::vector<std::vector<uint32_t>> total_fs;
1310 total_fs.reserve(M.elements.size() * 4);
1311 std::vector<std::tuple<uint32_t, uint32_t, uint32_t, uint32_t, uint32_t, uint32_t>> tempF;
1312 tempF.reserve(M.elements.size() * 4);
1313 std::vector<uint32_t> vs(3);
1314 for (const auto &ele : M_.elements)
1315 {
1316 for (int i = 0; i < 4; i++)
1317 {
1318 vs[0] = ele.vs[tet_faces[i][0]];
1319 vs[1] = ele.vs[tet_faces[i][1]];
1320 vs[2] = ele.vs[tet_faces[i][2]];
1321 total_fs.push_back(vs);
1322 std::sort(vs.begin(), vs.end());
1323 tempF.push_back(std::make_tuple(vs[0], vs[1], vs[2], ele.id * 4 + i, ele.id, i));
1324 }
1325 }
1326 std::sort(tempF.begin(), tempF.end());
1327 M_.faces.reserve(tempF.size() / 3);
1328 Face f;
1329 f.boundary = true;
1330 uint32_t F_num = 0;
1331 for (uint32_t i = 0; i < tempF.size(); ++i)
1332 {
1333 if (i == 0 || (i != 0 && (std::get<0>(tempF[i]) != std::get<0>(tempF[i - 1]) || std::get<1>(tempF[i]) != std::get<1>(tempF[i - 1]) || std::get<2>(tempF[i]) != std::get<2>(tempF[i - 1]))))
1334 {
1335 f.id = F_num;
1336 F_num++;
1337 f.vs = total_fs[std::get<3>(tempF[i])];
1338 M_.faces.push_back(f);
1339 }
1340 else if (i != 0 && (std::get<0>(tempF[i]) == std::get<0>(tempF[i - 1]) && std::get<1>(tempF[i]) == std::get<1>(tempF[i - 1]) && std::get<2>(tempF[i]) == std::get<2>(tempF[i - 1])))
1341 M_.faces[F_num - 1].boundary = false;
1342
1343 M_.elements[std::get<4>(tempF[i])].fs[std::get<5>(tempF[i])] = F_num - 1;
1344 }
1345
1346 build_connectivity(M_);
1347 orient_volume_mesh(M_);
1348 build_connectivity(M_);
1349
1350 M = M_;
1351
1352 double hmin = 10000, hmax = 0, havg = 0;
1353 for (const auto &e : M.edges)
1354 {
1355 Eigen::Vector3d v0 = M.points.col(e.vs[0]), v1 = M.points.col(e.vs[1]);
1356 double len = (v0 - v1).norm();
1357 if (len < hmin)
1358 hmin = len;
1359 if (len > hmax)
1360 hmax = len;
1361 havg += len;
1362 }
1363 havg /= M.edges.size();
1364 }
1365}
1366void MeshProcessing3D::straight_sweeping(const Mesh3DStorage &Mi, int sweep_coord, double height, int nlayer, Mesh3DStorage &Mo)
1367{
1368 if (sweep_coord > 2 || sweep_coord < 0)
1369 {
1370 logger().error("invalid sweeping direction!");
1371 return;
1372 }
1373 if (Mi.type != MeshType::H_SUR && Mi.type != MeshType::TRI && Mi.type != MeshType::QUA)
1374 {
1375 logger().error("invalid planar surface!");
1376 return;
1377 }
1378 if (height <= 0 || nlayer < 1)
1379 {
1380 logger().error("invalid height or number of layers!");
1381 return;
1382 }
1383
1384 Mo.points.resize(3, 0);
1385 Mo.vertices.clear();
1386 Mo.edges.clear();
1387 Mo.faces.clear();
1388 Mo.elements.clear();
1389 Mo.type = MeshType::HYB;
1390 // v, layers
1391 std::vector<std::vector<int>> Vlayers(nlayer + 1);
1392 std::vector<double> interval(3, 0);
1393 interval[sweep_coord] = height / nlayer;
1394
1395 for (int i = 0; i < nlayer + 1; i++)
1396 {
1397 std::vector<int> a_layer;
1398 for (auto v : Mi.vertices)
1399 {
1400 Vertex v_;
1401 v_.id = Mo.vertices.size();
1402 v_.v.resize(3);
1403 for (int j = 0; j < 3; j++)
1404 v_.v[j] = v.v[j] + i * interval[j];
1405
1406 Mo.vertices.push_back(v_);
1407 a_layer.push_back(v_.id);
1408 }
1409 Vlayers[i] = a_layer;
1410 }
1411 // f
1412 std::vector<std::vector<int>> Flayers(nlayer + 1);
1413 for (int i = 0; i < nlayer + 1; i++)
1414 {
1415 std::vector<int> a_layer;
1416 for (auto f : Mi.faces)
1417 {
1418 Face f_;
1419 f_.id = Mo.faces.size();
1420 for (auto vid : f.vs)
1421 f_.vs.push_back(Vlayers[i][vid]);
1422
1423 Mo.faces.push_back(f_);
1424 a_layer.push_back(f_.id);
1425 }
1426 Flayers[i] = a_layer;
1427 }
1428 // ef
1429 std::vector<std::vector<int>> EFlayers(nlayer);
1430 for (int i = 0; i < nlayer; i++)
1431 {
1432 std::vector<int> a_layer;
1433 for (auto e : Mi.edges)
1434 {
1435 Face f_;
1436 f_.id = Mo.faces.size();
1437 int v0 = e.vs[0], v1 = e.vs[1];
1438 f_.vs.push_back(Vlayers[i][v0]);
1439 f_.vs.push_back(Vlayers[i][v1]);
1440 f_.vs.push_back(Vlayers[i + 1][v1]);
1441 f_.vs.push_back(Vlayers[i + 1][v0]);
1442
1443 Mo.faces.push_back(f_);
1444 a_layer.push_back(f_.id);
1445 }
1446 EFlayers[i] = a_layer;
1447 }
1448 // ele
1449 for (int i = 0; i < nlayer; i++)
1450 {
1451 for (auto f : Mi.faces)
1452 {
1453 Element ele;
1454 ele.id = Mo.elements.size();
1455 ele.hex = (f.vs.size() == 4);
1456
1457 ele.fs.push_back(Flayers[i][f.id]);
1458 ele.fs.push_back(Flayers[i + 1][f.id]);
1459 for (auto eid : f.es)
1460 ele.fs.push_back(EFlayers[i][eid]);
1461
1462 ele.fs_flag.resize(ele.fs.size(), false);
1463
1464 ele.v_in_Kernel.resize(3, 0);
1465 int nv = 0;
1466 for (int j = 0; j < 2; j++)
1467 {
1468 nv += Mo.faces[ele.fs[j]].vs.size();
1469 for (auto vid : Mo.faces[ele.fs[j]].vs)
1470 for (int k = 0; k < 3; k++)
1471 ele.v_in_Kernel[k] += Mo.vertices[vid].v[k];
1472 }
1473 for (int k = 0; k < 3; k++)
1474 ele.v_in_Kernel[k] /= nv;
1475
1476 Mo.elements.push_back(ele);
1477 }
1478 }
1479
1480 Mo.points.resize(3, Mo.vertices.size());
1481 for (auto v : Mo.vertices)
1482 for (int j = 0; j < 3; j++)
1483 Mo.points(j, v.id) = v.v[j];
1484
1485 build_connectivity(Mo);
1486 orient_volume_mesh(Mo);
1487 build_connectivity(Mo);
1488}
1489
1490void MeshProcessing3D::orient_surface_mesh(Mesh3DStorage &hmi)
1491{
1492
1493 vector<bool> flag(hmi.faces.size(), true);
1494 flag[0] = false;
1495
1496 std::queue<uint32_t> pf_temp;
1497 pf_temp.push(0);
1498 while (!pf_temp.empty())
1499 {
1500 uint32_t fid = pf_temp.front();
1501 pf_temp.pop();
1502 for (auto eid : hmi.faces[fid].es)
1503 for (auto nfid : hmi.edges[eid].neighbor_fs)
1504 {
1505 if (!flag[nfid])
1506 continue;
1507 uint32_t v0 = hmi.edges[eid].vs[0], v1 = hmi.edges[eid].vs[1];
1508 int32_t v0_pos = std::find(hmi.faces[fid].vs.begin(), hmi.faces[fid].vs.end(), v0) - hmi.faces[fid].vs.begin();
1509 int32_t v1_pos = std::find(hmi.faces[fid].vs.begin(), hmi.faces[fid].vs.end(), v1) - hmi.faces[fid].vs.begin();
1510
1511 if ((v0_pos + 1) % hmi.faces[fid].vs.size() != v1_pos)
1512 swap(v0, v1);
1513
1514 int32_t v0_pos_ = std::find(hmi.faces[nfid].vs.begin(), hmi.faces[nfid].vs.end(), v0) - hmi.faces[nfid].vs.begin();
1515 int32_t v1_pos_ = std::find(hmi.faces[nfid].vs.begin(), hmi.faces[nfid].vs.end(), v1) - hmi.faces[nfid].vs.begin();
1516
1517 if ((v0_pos_ + 1) % hmi.faces[nfid].vs.size() == v1_pos_)
1518 std::reverse(hmi.faces[nfid].vs.begin(), hmi.faces[nfid].vs.end());
1519
1520 pf_temp.push(nfid);
1521 flag[nfid] = false;
1522 }
1523 }
1524 double res = 0;
1525 Vector3d ori;
1526 ori.setZero();
1527 for (auto f : hmi.faces)
1528 {
1529 auto &fvs = f.vs;
1530 Vector3d center;
1531 center.setZero();
1532 for (auto vid : fvs)
1533 center += hmi.points.col(vid);
1534 center /= fvs.size();
1535
1536 for (uint32_t j = 0; j < fvs.size(); j++)
1537 {
1538 Vector3d x = hmi.points.col(fvs[j]) - ori, y = hmi.points.col(fvs[(j + 1) % fvs.size()]) - ori, z = center - ori;
1539 res += -((x[0] * y[1] * z[2] + x[1] * y[2] * z[0] + x[2] * y[0] * z[1]) - (x[2] * y[1] * z[0] + x[1] * y[0] * z[2] + x[0] * y[2] * z[1]));
1540 }
1541 }
1542 if (res > 0)
1543 {
1544 for (uint32_t i = 0; i < hmi.faces.size(); i++)
1545 std::reverse(hmi.faces[i].vs.begin(), hmi.faces[i].vs.end());
1546 }
1547}
1548void MeshProcessing3D::orient_volume_mesh(Mesh3DStorage &hmi)
1549{
1550 // surface orienting
1551 Mesh3DStorage M_sur;
1552 M_sur.type = MeshType::H_SUR;
1553 int bvn = 0;
1554 for (auto v : hmi.vertices)
1555 if (v.boundary)
1556 bvn++;
1557 M_sur.points.resize(3, bvn);
1558 bvn = 0;
1559 vector<int> V_map(hmi.vertices.size(), -1), V_map_reverse;
1560 for (auto v : hmi.vertices)
1561 if (v.boundary)
1562 {
1563 M_sur.points.col(bvn++) = hmi.points.col(v.id);
1564 Vertex v_;
1565 v_.id = M_sur.vertices.size();
1566 M_sur.vertices.push_back(v_);
1567
1568 V_map[v.id] = v_.id;
1569 V_map_reverse.push_back(v.id);
1570 }
1571 for (auto f : hmi.faces)
1572 if (f.boundary)
1573 {
1574 for (int j = 0; j < f.vs.size(); j++)
1575 {
1576 f.id = M_sur.faces.size();
1577 f.es.clear();
1578 f.neighbor_hs.clear();
1579 f.vs[j] = V_map[f.vs[j]];
1580 M_sur.vertices[f.vs[j]].neighbor_fs.push_back(f.id);
1581 }
1582 M_sur.faces.push_back(f);
1583 }
1584 build_connectivity(M_sur);
1585 orient_surface_mesh(M_sur);
1586
1587 int fn_ = 0;
1588 for (auto &f : hmi.faces)
1589 if (f.boundary)
1590 {
1591 for (int j = 0; j < f.vs.size(); j++)
1592 f.vs[j] = V_map_reverse[M_sur.faces[fn_].vs[j]];
1593 fn_++;
1594 }
1595 // volume orienting
1596 vector<bool> F_tag(hmi.faces.size(), true);
1597 std::vector<short> F_visit(hmi.faces.size(), 0); // 0 un-visited, 1 visited once, 2 visited twice
1598 for (uint32_t j = 0; j < hmi.faces.size(); j++)
1599 if (hmi.faces[j].boundary)
1600 {
1601 F_visit[j]++;
1602 }
1603 std::vector<bool> F_state(hmi.faces.size(), false); // false is the reverse direction, true is the same direction
1604 std::vector<bool> P_visit(hmi.elements.size(), false);
1605 while (true)
1606 {
1607 std::vector<uint32_t> candidates;
1608 for (uint32_t j = 0; j < F_visit.size(); j++)
1609 if (F_visit[j] == 1)
1610 candidates.push_back(j);
1611 if (!candidates.size())
1612 break;
1613 for (auto ca : candidates)
1614 {
1615 if (F_visit[ca] == 2)
1616 continue;
1617 uint32_t pid = hmi.faces[ca].neighbor_hs[0];
1618 if (P_visit[pid])
1619 if (hmi.faces[ca].neighbor_hs.size() == 2)
1620 pid = hmi.faces[ca].neighbor_hs[1];
1621 if (P_visit[pid])
1622 {
1623 logger().error("bug");
1624 }
1625 auto &fs = hmi.elements[pid].fs;
1626 for (auto fid : fs)
1627 F_tag[fid] = false;
1628
1629 uint32_t start_f = ca;
1630 F_tag[start_f] = true;
1631 F_visit[ca]++;
1632 if (F_state[ca])
1633 F_state[ca] = false;
1634 else
1635 F_state[ca] = true;
1636
1637 std::queue<uint32_t> pf_temp;
1638 pf_temp.push(start_f);
1639 while (!pf_temp.empty())
1640 {
1641 uint32_t fid = pf_temp.front();
1642 pf_temp.pop();
1643 for (auto eid : hmi.faces[fid].es)
1644 for (auto nfid : hmi.edges[eid].neighbor_fs)
1645 {
1646
1647 if (F_tag[nfid])
1648 continue;
1649 uint32_t v0 = hmi.edges[eid].vs[0], v1 = hmi.edges[eid].vs[1];
1650 int32_t v0_pos = std::find(hmi.faces[fid].vs.begin(), hmi.faces[fid].vs.end(), v0) - hmi.faces[fid].vs.begin();
1651 int32_t v1_pos = std::find(hmi.faces[fid].vs.begin(), hmi.faces[fid].vs.end(), v1) - hmi.faces[fid].vs.begin();
1652
1653 if ((v0_pos + 1) % hmi.faces[fid].vs.size() != v1_pos)
1654 std::swap(v0, v1);
1655
1656 int32_t v0_pos_ = std::find(hmi.faces[nfid].vs.begin(), hmi.faces[nfid].vs.end(), v0) - hmi.faces[nfid].vs.begin();
1657 int32_t v1_pos_ = std::find(hmi.faces[nfid].vs.begin(), hmi.faces[nfid].vs.end(), v1) - hmi.faces[nfid].vs.begin();
1658
1659 if (F_state[fid])
1660 {
1661 if ((v0_pos_ + 1) % hmi.faces[nfid].vs.size() == v1_pos_)
1662 F_state[nfid] = false;
1663 else
1664 F_state[nfid] = true;
1665 }
1666 else if (!F_state[fid])
1667 {
1668 if ((v0_pos_ + 1) % hmi.faces[nfid].vs.size() == v1_pos_)
1669 F_state[nfid] = true;
1670 else
1671 F_state[nfid] = false;
1672 }
1673
1674 F_visit[nfid]++;
1675
1676 pf_temp.push(nfid);
1677 F_tag[nfid] = true;
1678 }
1679 }
1680 P_visit[pid] = true;
1681 for (uint32_t j = 0; j < fs.size(); j++)
1682 hmi.elements[pid].fs_flag[j] = F_state[fs[j]];
1683 }
1684 }
1685}
1686void MeshProcessing3D::ele_subdivison_levels(const Mesh3DStorage &hmi, std::vector<int> &Ls)
1687{
1688
1689 Ls.clear();
1690 Ls.resize(hmi.elements.size(), 1);
1691 std::vector<double> volumes(hmi.elements.size(), 0);
1692
1693 auto compute_volume = [&](const int id, double &vol) {
1694 Vector3d ori;
1695 ori.setZero();
1696 for (auto f : hmi.elements[id].fs)
1697 {
1698 auto &fvs = hmi.faces[f].vs;
1699 Vector3d center;
1700 center.setZero();
1701 for (auto vid : fvs)
1702 center += hmi.points.col(vid);
1703 center /= fvs.size();
1704
1705 for (uint32_t j = 0; j < fvs.size(); j++)
1706 {
1707 Vector3d x = hmi.points.col(fvs[j]) - ori, y = hmi.points.col(fvs[(j + 1) % fvs.size()]) - ori, z = center - ori;
1708 vol += -((x[0] * y[1] * z[2] + x[1] * y[2] * z[0] + x[2] * y[0] * z[1]) - (x[2] * y[1] * z[0] + x[1] * y[0] * z[2] + x[0] * y[2] * z[1]));
1709 }
1710 }
1711 vol = std::abs(vol);
1712 };
1713
1714 for (const auto &ele : hmi.elements)
1715 compute_volume(ele.id, volumes[ele.id]);
1716
1717 double ave_volume = 0;
1718 for (const auto &v : volumes)
1719 ave_volume += v;
1720 ave_volume /= volumes.size();
1721 for (int i = 0; i < Ls.size(); i++)
1722 if (!hmi.elements[i].hex)
1723 {
1724 Ls[i] = volumes[i] / ave_volume;
1725 if (Ls[i] < 1)
1726 Ls[i] = 1;
1727 }
1728}
1729
1730// template<typename T>
1731// void MeshProcessing3D::set_intersection_own(const std::vector<T> &A, const std::vector<T> &B, std::vector<T> &C, const int &num){
1732void MeshProcessing3D::set_intersection_own(const std::vector<uint32_t> &A, const std::vector<uint32_t> &B, std::array<uint32_t, 2> &C, int &num)
1733{
1734 // void MeshProcessing3D::set_intersection_own( std::vector<uint32_t> &A, std::vector<uint32_t> &B, std::vector<uint32_t> &C, int &num)
1735 // C.resize(num);
1736 int n = 0;
1737 for (auto &a : A)
1738 {
1739 for (auto &b : B)
1740 {
1741 if (a == b)
1742 {
1743 C[n++] = a;
1744 if (n == num)
1745 break;
1746 }
1747 }
1748 if (n == num)
1749 break;
1750 }
1751}
Jacobian_Precision
#define Jacobian_Precision
Definition MeshProcessing3D.hpp:14
y
int y
Definition SplineBasis3d.cpp:55
z
int z
Definition SplineBasis3d.cpp:55
x
int x
Definition SplineBasis3d.cpp:55
polyfem::mesh::Mesh3DStorage::elements
std::vector< Element > elements
Definition Mesh3DStorage.hpp:75
polyfem::mesh::Mesh3DStorage::vertices
std::vector< Vertex > vertices
Definition Mesh3DStorage.hpp:72
polyfem::mesh::MeshProcessing3D::reorder_hex_mesh_propogation
void reorder_hex_mesh_propogation(Mesh3DStorage &hmi)
Definition MeshProcessing3D.cpp:493
polyfem::mesh::MeshProcessing3D::global_orientation_hexes
void global_orientation_hexes(Mesh3DStorage &hmi)
Definition MeshProcessing3D.cpp:676
polyfem::mesh::MeshProcessing3D::a_jacobian
double a_jacobian(Eigen::Vector3d &v0, Eigen::Vector3d &v1, Eigen::Vector3d &v2, Eigen::Vector3d &v3)
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:1732
polyfem::mesh::MeshProcessing3D::refine_red_refinement_tet
void refine_red_refinement_tet(Mesh3DStorage &M, int iter)
Definition MeshProcessing3D.cpp:1117
polyfem::mesh::MeshProcessing3D::straight_sweeping
void straight_sweeping(const Mesh3DStorage &Mi, int sweep_coord, double height, int nlayer, Mesh3DStorage &Mo)
Definition MeshProcessing3D.cpp:1366
polyfem::mesh::MeshProcessing3D::orient_surface_mesh
void orient_surface_mesh(Mesh3DStorage &hmi)
Definition MeshProcessing3D.cpp:1490
polyfem::mesh::MeshProcessing3D::build_connectivity
void build_connectivity(Mesh3DStorage &hmi)
Definition MeshProcessing3D.cpp:19
polyfem::mesh::MeshProcessing3D::orient_volume_mesh
void orient_volume_mesh(Mesh3DStorage &hmi)
Definition MeshProcessing3D.cpp:1548
polyfem::mesh::MeshProcessing3D::ele_subdivison_levels
void ele_subdivison_levels(const Mesh3DStorage &hmi, std::vector< int > &Ls)
Definition MeshProcessing3D.cpp:1686
polyfem::mesh::MeshProcessing3D::scaled_jacobian
bool scaled_jacobian(Mesh3DStorage &hmi, Mesh_Quality &mq)
Definition MeshProcessing3D.cpp:606
polyfem::mesh::MeshProcessing3D::refine_catmul_clark_polar
void refine_catmul_clark_polar(Mesh3DStorage &M, int iter, bool reverse, std::vector< int > &Parents)
Definition MeshProcessing3D.cpp:711
polyfem::mesh::find
int find(const Eigen::VectorXi &vec, int x)
Definition NCMesh2D.cpp:341
polyfem::mesh::signed_volume
double signed_volume(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F)
Compute the signed volume of a surface mesh.
Definition MeshUtils.cpp:677
polyfem::logger
spdlog::logger & logger()
Retrieves the current logger.
Definition Logger.cpp:44
polyfem::mesh::Element::v_in_Kernel
std::vector< double > v_in_Kernel
Definition Mesh3DStorage.hpp:54
polyfem::mesh::Element::fs
std::vector< uint32_t > fs
Definition Mesh3DStorage.hpp:51
polyfem::mesh::Element::vs
std::vector< uint32_t > vs
Definition Mesh3DStorage.hpp:49
polyfem::mesh::Element::fs_flag
std::vector< bool > fs_flag
Definition Mesh3DStorage.hpp:52
polyfem::mesh::Face::vs
std::vector< uint32_t > vs
Definition Mesh3DStorage.hpp:36
polyfem::mesh::Mesh_Quality
Definition Mesh3DStorage.hpp:205
polyfem::mesh::Mesh_Quality::deviation_Jacobian
double deviation_Jacobian
Definition Mesh3DStorage.hpp:209
polyfem::mesh::Mesh_Quality::min_Jacobian
double min_Jacobian
Definition Mesh3DStorage.hpp:207
polyfem::mesh::Mesh_Quality::H_Js
Eigen::VectorXd H_Js
Definition Mesh3DStorage.hpp:211
polyfem::mesh::Mesh_Quality::V_Js
Eigen::VectorXd V_Js
Definition Mesh3DStorage.hpp:210
polyfem::mesh::Mesh_Quality::ave_Jacobian
double ave_Jacobian
Definition Mesh3DStorage.hpp:208
polyfem::mesh::Vertex::v
std::vector< double > v
Definition Mesh3DStorage.hpp:14