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