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
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 // cout << "correct orientation " << endl;
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 // cout << "m1 jacobian " << mq1.min_Jacobian << " " << mq1.ave_Jacobian << endl;
591 if (mq1.min_Jacobian > 0)
592 continue;
593 m2.elements = m1.elements;
594 for (auto &h : m2.elements)
595 {
596 swap(h.vs[1], h.vs[3]);
597 swap(h.vs[5], h.vs[7]);
598 }
599 scaled_jacobian(m2, mq2);
600 // cout << "m2 jacobian " << mq2.min_Jacobian << " " << mq2.ave_Jacobian << endl;
601 if (mq2.ave_Jacobian > mq1.ave_Jacobian)
602 {
603 for (auto &h : m2.elements)
604 hmi.elements[h.id] = h;
605 }
606 }
607}
608bool MeshProcessing3D::scaled_jacobian(Mesh3DStorage &hmi, Mesh_Quality &mq)
609{
610 if (hmi.type != MeshType::HEX)
611 return false;
612
613 mq.ave_Jacobian = 0;
614 mq.min_Jacobian = 1;
615 mq.deviation_Jacobian = 0;
616 mq.V_Js.resize(hmi.elements.size() * 8);
617 mq.V_Js.setZero();
618 mq.H_Js.resize(hmi.elements.size());
619 mq.H_Js.setZero();
620
621 for (uint32_t i = 0; i < hmi.elements.size(); i++)
622 {
623 double hex_minJ = 1;
624 for (uint32_t j = 0; j < 8; j++)
625 {
626 uint32_t v0, v1, v2, v3;
627 v0 = hex_tetra_table[j][0];
628 v1 = hex_tetra_table[j][1];
629 v2 = hex_tetra_table[j][2];
630 v3 = hex_tetra_table[j][3];
631
632 Vector3d c0 = hmi.points.col(hmi.elements[i].vs[v0]);
633 Vector3d c1 = hmi.points.col(hmi.elements[i].vs[v1]);
634 Vector3d c2 = hmi.points.col(hmi.elements[i].vs[v2]);
635 Vector3d c3 = hmi.points.col(hmi.elements[i].vs[v3]);
636
637 double jacobian_value = a_jacobian(c0, c1, c2, c3);
638
639 if (hex_minJ > jacobian_value)
640 hex_minJ = jacobian_value;
641
642 uint32_t id = 8 * i + j;
643 mq.V_Js[id] = jacobian_value;
644 }
645 mq.H_Js[i] = hex_minJ;
646 mq.ave_Jacobian += hex_minJ;
647 if (mq.min_Jacobian > hex_minJ)
648 mq.min_Jacobian = hex_minJ;
649 }
650 mq.ave_Jacobian /= hmi.elements.size();
651 for (int i = 0; i < mq.H_Js.size(); i++)
652 mq.deviation_Jacobian += (mq.H_Js[i] - mq.ave_Jacobian) * (mq.H_Js[i] - mq.ave_Jacobian);
653 mq.deviation_Jacobian /= hmi.elements.size();
654
655 return true;
656}
657double MeshProcessing3D::a_jacobian(Vector3d &v0, Vector3d &v1, Vector3d &v2, Vector3d &v3)
658{
659 Matrix3d Jacobian;
660
661 Jacobian.col(0) = v1 - v0;
662 Jacobian.col(1) = v2 - v0;
663 Jacobian.col(2) = v3 - v0;
664
665 double norm1 = Jacobian.col(0).norm();
666 double norm2 = Jacobian.col(1).norm();
667 double norm3 = Jacobian.col(2).norm();
668
669 double scaled_jacobian = Jacobian.determinant();
670 if (std::abs(norm1) < Jacobian_Precision || std::abs(norm2) < Jacobian_Precision || std::abs(norm3) < Jacobian_Precision)
671 {
672 // std::cout << "Potential Bug, check!" << endl; //system("PAUSE");
673 return scaled_jacobian;
674 }
675 scaled_jacobian /= norm1 * norm2 * norm3;
676 return scaled_jacobian;
677}
678
679void MeshProcessing3D::global_orientation_hexes(Mesh3DStorage &hmi)
680{
681 Mesh3DStorage mesh;
682 mesh.type = MeshType::HEX;
683 mesh.points = hmi.points;
684 mesh.vertices.resize(hmi.vertices.size());
685 for (const auto &v : hmi.vertices)
686 {
687 Vertex v_;
688 v_.id = v.id;
689 v_.v = v.v;
690 mesh.vertices[v.id] = v;
691 }
692 vector<int> Ele_map(hmi.elements.size(), -1), Ele_map_reverse;
693 for (const auto &ele : hmi.elements)
694 {
695 if (!ele.hex)
696 continue;
697 Element ele_;
698 ele_.id = mesh.elements.size();
699 ele_.vs = ele.vs;
700 for (auto vid : ele_.vs)
701 mesh.vertices[vid].neighbor_hs.push_back(ele_.id);
702 mesh.elements.push_back(ele_);
703
704 Ele_map[ele.id] = ele_.id;
705 Ele_map_reverse.push_back(ele.id);
706 }
707
708 build_connectivity(mesh);
709 reorder_hex_mesh_propogation(mesh);
710
711 for (const auto &ele : mesh.elements)
712 hmi.elements[Ele_map_reverse[ele.id]].vs = ele.vs;
713}
714void MeshProcessing3D::refine_catmul_clark_polar(Mesh3DStorage &M, int iter, bool reverse, std::vector<int> &Parents)
715{
716
717 for (int i = 0; i < iter; i++)
718 {
719
720 std::vector<int> Refinement_Levels;
721 ele_subdivison_levels(M, Refinement_Levels);
722
723 Mesh3DStorage M_;
724 M_.type = MeshType::HYB;
725
726 vector<int> E2V(M.edges.size()), F2V(M.faces.size()), Ele2V(M.elements.size());
727
728 int vn = 0;
729 for (auto v : M.vertices)
730 {
731 Vertex v_;
732 v_.id = vn++;
733 v_.v.resize(3);
734 for (int j = 0; j < 3; j++)
735 v_.v[j] = M.points(j, v.id);
736 M_.vertices.push_back(v_);
737 }
738
739 for (auto e : M.edges)
740 {
741 Vertex v;
742 v.id = vn++;
743 v.v.resize(3);
744
745 Vector3d center;
746 center.setZero();
747 for (auto vid : e.vs)
748 center += M.points.col(vid);
749 center /= e.vs.size();
750 for (int j = 0; j < 3; j++)
751 v.v[j] = center[j];
752
753 M_.vertices.push_back(v);
754 E2V[e.id] = v.id;
755 }
756 for (auto f : M.faces)
757 {
758 Vertex v;
759 v.id = vn++;
760 v.v.resize(3);
761
762 Vector3d center;
763 center.setZero();
764 for (auto vid : f.vs)
765 center += M.points.col(vid);
766 center /= f.vs.size();
767 for (int j = 0; j < 3; j++)
768 v.v[j] = center[j];
769
770 M_.vertices.push_back(v);
771 F2V[f.id] = v.id;
772 }
773 for (auto ele : M.elements)
774 {
775 if (!ele.hex)
776 continue;
777 Vertex v;
778 v.id = vn++;
779 v.v = ele.v_in_Kernel;
780
781 M_.vertices.push_back(v);
782 Ele2V[ele.id] = v.id;
783 }
784 // new elements
785 std::vector<std::vector<uint32_t>> total_fs;
786 total_fs.reserve(M.elements.size() * 8 * 6);
787 std::vector<std::tuple<uint32_t, uint32_t, uint32_t, uint32_t, uint32_t, uint32_t, uint32_t>> tempF;
788 tempF.reserve(M.elements.size() * 8 * 6);
789 std::vector<uint32_t> vs(4);
790
791 int elen = 0, fn = 0;
792 for (auto ele : M.elements)
793 {
794 if (ele.hex)
795 {
796 for (auto vid : ele.vs)
797 {
798 // top 4 vs
799 vector<int> top_vs(4);
800 top_vs[0] = vid;
801 int fid = -1;
802 for (auto nfid : M.vertices[vid].neighbor_fs)
803 if (find(ele.fs.begin(), ele.fs.end(), nfid) != ele.fs.end())
804 {
805 fid = nfid;
806 break;
807 }
808 assert(fid != -1);
809 top_vs[2] = F2V[fid];
810
811 int v_ind = find(M.faces[fid].vs.begin(), M.faces[fid].vs.end(), vid) - M.faces[fid].vs.begin();
812 int e_pre = M.faces[fid].es[(v_ind - 1 + 4) % 4];
813 int e_aft = M.faces[fid].es[v_ind];
814 top_vs[1] = E2V[e_pre];
815 top_vs[3] = E2V[e_aft];
816 // bottom 4 vs
817 vector<int> bottom_vs(4);
818
819 auto nvs = M.vertices[vid].neighbor_vs;
820 int e_per = -1;
821 for (auto nvid : nvs)
822 if (find(ele.vs.begin(), ele.vs.end(), nvid) != ele.vs.end())
823 {
824 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])
825 {
826 vector<uint32_t> sharedes, es0 = M.vertices[vid].neighbor_es, es1 = M.vertices[nvid].neighbor_es;
827 sort(es0.begin(), es0.end());
828 sort(es1.begin(), es1.end());
829 set_intersection(es0.begin(), es0.end(), es1.begin(), es1.end(), back_inserter(sharedes));
830 assert(sharedes.size());
831 e_per = sharedes[0];
832 break;
833 }
834 }
835
836 assert(e_per != -1);
837 bottom_vs[0] = E2V[e_per];
838 bottom_vs[2] = Ele2V[ele.id];
839
840 int f_pre = -1;
841 vector<uint32_t> sharedfs, fs0 = M.edges[e_pre].neighbor_fs, fs1 = M.edges[e_per].neighbor_fs;
842 sort(fs0.begin(), fs0.end());
843 sort(fs1.begin(), fs1.end());
844 set_intersection(fs0.begin(), fs0.end(), fs1.begin(), fs1.end(), back_inserter(sharedfs));
845 for (auto sfid : sharedfs)
846 if (find(ele.fs.begin(), ele.fs.end(), sfid) != ele.fs.end())
847 {
848 f_pre = sfid;
849 break;
850 }
851 assert(f_pre != -1);
852
853 int f_aft = -1;
854 sharedfs.clear();
855 fs0 = M.edges[e_aft].neighbor_fs;
856 fs1 = M.edges[e_per].neighbor_fs;
857 sort(fs0.begin(), fs0.end());
858 sort(fs1.begin(), fs1.end());
859 set_intersection(fs0.begin(), fs0.end(), fs1.begin(), fs1.end(), back_inserter(sharedfs));
860 for (auto sfid : sharedfs)
861 if (find(ele.fs.begin(), ele.fs.end(), sfid) != ele.fs.end())
862 {
863 f_aft = sfid;
864 break;
865 }
866 assert(f_aft != -1);
867
868 bottom_vs[1] = F2V[f_pre];
869 bottom_vs[3] = F2V[f_aft];
870
871 vector<int> ele_vs = top_vs;
872 ele_vs.insert(ele_vs.end(), bottom_vs.begin(), bottom_vs.end());
873 // fs
874 for (short j = 0; j < 6; j++)
875 {
876 for (short k = 0; k < 4; k++)
877 vs[k] = ele_vs[hex_face_table[j][k]];
878 total_fs.push_back(vs);
879 std::sort(vs.begin(), vs.end());
880 tempF.push_back(std::make_tuple(vs[0], vs[1], vs[2], vs[3], fn++, elen, j));
881 }
882 // new ele
883 Element ele_;
884 ele_.id = elen++;
885 ele_.fs.resize(6, -1);
886 ele_.fs_flag.resize(6, 1);
887
888 ele_.hex = true;
889 ele_.v_in_Kernel.resize(3);
890
891 Vector3d center;
892 center.setZero();
893 for (auto evid : ele_vs)
894 for (int j = 0; j < 3; j++)
895 center[j] += M_.vertices[evid].v[j];
896 center /= ele_vs.size();
897 for (int j = 0; j < 3; j++)
898 ele_.v_in_Kernel[j] = center[j];
899
900 M_.elements.push_back(ele_);
901 Parents.push_back(ele.id);
902 }
903 }
904 else
905 {
906 int level = Refinement_Levels[ele.id];
907 if (reverse)
908 level = 1;
909 // local_V2V
910 std::vector<std::vector<int>> local_V2Vs;
911 std::map<int, int> local_vi_map;
912 for (auto vid : ele.vs)
913 {
914 std::vector<int> v2v;
915 v2v.push_back(vid);
916 for (int r = 0; r < level; r++)
917 {
918 Vertex v_;
919 v_.id = vn++;
920 v_.v.resize(3);
921
922 for (int j = 0; j < 3; j++)
923 v_.v[j] = M_.vertices[vid].v[j] + (ele.v_in_Kernel[j] - M_.vertices[vid].v[j]) * (r + 1.0) / (double)(level + 1);
924 // cout << "before: "<<v_.v[0] << " " << v_.v[1] << " " << v_.v[2] << endl;
925 if (reverse)
926 {
927 for (int j = 0; j < 3; j++)
928 v_.v[j] = M_.vertices[vid].v[j] + (M_.vertices[vid].v[j] - ele.v_in_Kernel[j]) * (r + 1.0) / (double)(level + 1);
929 // cout << "after: " << v_.v[0] << " " << v_.v[1] << " " << v_.v[2] << endl;
930 }
931 M_.vertices.push_back(v_);
932 v2v.push_back(v_.id);
933 }
934 local_vi_map[vid] = local_V2Vs.size();
935 local_V2Vs.push_back(v2v);
936 }
937 // local_E2V
938 vector<uint32_t> es;
939 for (auto fid : ele.fs)
940 es.insert(es.end(), M.faces[fid].es.begin(), M.faces[fid].es.end());
941 sort(es.begin(), es.end());
942 es.erase(unique(es.begin(), es.end()), es.end());
943
944 std::vector<std::vector<int>> local_E2Vs;
945 std::map<int, int> local_ei_map;
946 for (auto eid : es)
947 {
948 std::vector<int> e2v;
949 e2v.push_back(E2V[eid]);
950 for (int r = 0; r < level; r++)
951 {
952 Vertex v_;
953 v_.id = vn++;
954 v_.v.resize(3);
955
956 Vector3d center;
957 center.setZero();
958 for (auto vid : M.edges[eid].vs)
959 for (int j = 0; j < 3; j++)
960 center[j] += M_.vertices[local_V2Vs[local_vi_map[vid]][r + 1]].v[j];
961 center /= M.edges[eid].vs.size();
962 for (int j = 0; j < 3; j++)
963 v_.v[j] = center[j];
964
965 M_.vertices.push_back(v_);
966 e2v.push_back(v_.id);
967 }
968 local_ei_map[eid] = local_E2Vs.size();
969 local_E2Vs.push_back(e2v);
970 }
971 // local_F2V
972 std::vector<std::vector<int>> local_F2Vs;
973 std::map<int, int> local_fi_map;
974 for (auto fid : ele.fs)
975 {
976 std::vector<int> f2v;
977 f2v.push_back(F2V[fid]);
978 for (int r = 0; r < level; r++)
979 {
980 Vertex v_;
981 v_.id = vn++;
982 v_.v.resize(3);
983
984 Vector3d center;
985 center.setZero();
986 for (auto vid : M.faces[fid].vs)
987 for (int j = 0; j < 3; j++)
988 center[j] += M_.vertices[local_V2Vs[local_vi_map[vid]][r + 1]].v[j];
989 center /= M.faces[fid].vs.size();
990 for (int j = 0; j < 3; j++)
991 v_.v[j] = center[j];
992
993 M_.vertices.push_back(v_);
994 f2v.push_back(v_.id);
995 }
996 local_fi_map[fid] = local_F2Vs.size();
997 local_F2Vs.push_back(f2v);
998 }
999 // polyhedron fs
1000 int local_fn = 0;
1001 for (auto fid : ele.fs)
1002 {
1003 auto &fvs = M.faces[fid].vs;
1004 auto &fes = M.faces[fid].es;
1005 int fvn = M.faces[fid].vs.size();
1006 for (uint32_t j = 0; j < fvn; j++)
1007 {
1008 vs[0] = local_E2Vs[local_ei_map[fes[(j - 1 + fvn) % fvn]]][level];
1009 vs[1] = local_V2Vs[local_vi_map[fvs[j]]][level];
1010 vs[2] = local_E2Vs[local_ei_map[fes[j]]][level];
1011 vs[3] = local_F2Vs[local_fi_map[fid]][level];
1012
1013 if (reverse)
1014 {
1015 vs[0] = local_E2Vs[local_ei_map[fes[(j - 1 + fvn) % fvn]]][0];
1016 vs[1] = local_V2Vs[local_vi_map[fvs[j]]][0];
1017 vs[2] = local_E2Vs[local_ei_map[fes[j]]][0];
1018 vs[3] = local_F2Vs[local_fi_map[fid]][0];
1019 }
1020
1021 total_fs.push_back(vs);
1022 std::sort(vs.begin(), vs.end());
1023 tempF.push_back(std::make_tuple(vs[0], vs[1], vs[2], vs[3], fn++, elen, local_fn++));
1024 }
1025 }
1026 // polyhedron
1027 Element ele_;
1028 ele_.id = elen++;
1029 ele_.fs.resize(local_fn, -1);
1030 ele_.fs_flag.resize(local_fn, 1);
1031
1032 ele_.hex = false;
1033 ele_.v_in_Kernel = ele.v_in_Kernel;
1034 M_.elements.push_back(ele_);
1035 Parents.push_back(ele.id);
1036 // hex
1037 for (int r = 0; r < level; r++)
1038 {
1039 for (auto fid : ele.fs)
1040 {
1041 auto &fvs = M.faces[fid].vs;
1042 auto &fes = M.faces[fid].es;
1043 int fvn = M.faces[fid].vs.size();
1044 for (uint32_t j = 0; j < fvn; j++)
1045 {
1046 vector<int> ele_vs(8);
1047 ele_vs[0] = local_E2Vs[local_ei_map[fes[(j - 1 + fvn) % fvn]]][r + 1];
1048 ele_vs[1] = local_V2Vs[local_vi_map[fvs[j]]][r + 1];
1049 ele_vs[2] = local_E2Vs[local_ei_map[fes[j]]][r + 1];
1050 ele_vs[3] = local_F2Vs[local_fi_map[fid]][r + 1];
1051
1052 ele_vs[4] = local_E2Vs[local_ei_map[fes[(j - 1 + fvn) % fvn]]][r];
1053 ele_vs[5] = local_V2Vs[local_vi_map[fvs[j]]][r];
1054 ele_vs[6] = local_E2Vs[local_ei_map[fes[j]]][r];
1055 ele_vs[7] = local_F2Vs[local_fi_map[fid]][r];
1056
1057 // fs
1058 for (short j = 0; j < 6; j++)
1059 {
1060 for (short k = 0; k < 4; k++)
1061 vs[k] = ele_vs[hex_face_table[j][k]];
1062 total_fs.push_back(vs);
1063 std::sort(vs.begin(), vs.end());
1064 tempF.push_back(std::make_tuple(vs[0], vs[1], vs[2], vs[3], fn++, elen, j));
1065 }
1066 // hex
1067 Element ele_;
1068 ele_.id = elen++;
1069 ele_.fs.resize(6, -1);
1070 ele_.fs_flag.resize(6, 1);
1071 ele_.hex = true;
1072 ele_.v_in_Kernel.resize(3);
1073
1074 Vector3d center;
1075 center.setZero();
1076 for (auto vid : ele_vs)
1077 for (int j = 0; j < 3; j++)
1078 center[j] += M_.vertices[vid].v[j];
1079 center /= ele_vs.size();
1080 for (int j = 0; j < 3; j++)
1081 ele_.v_in_Kernel[j] = center[j];
1082 M_.elements.push_back(ele_);
1083 Parents.push_back(ele.id);
1084 }
1085 }
1086 }
1087 }
1088 }
1089 // Fs
1090 std::sort(tempF.begin(), tempF.end());
1091 M_.faces.reserve(tempF.size() / 3);
1092 Face f;
1093 f.boundary = true;
1094 uint32_t F_num = 0;
1095 for (uint32_t i = 0; i < tempF.size(); ++i)
1096 {
1097 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]))))
1098 {
1099 f.id = F_num;
1100 F_num++;
1101 f.vs = total_fs[std::get<4>(tempF[i])];
1102 M_.faces.push_back(f);
1103 }
1104 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])))
1105 M_.faces[F_num - 1].boundary = false;
1106
1107 M_.elements[std::get<5>(tempF[i])].fs[std::get<6>(tempF[i])] = F_num - 1;
1108 }
1109
1110 M_.points.resize(3, M_.vertices.size());
1111 for (auto v : M_.vertices)
1112 for (int j = 0; j < 3; j++)
1113 M_.points(j, v.id) = v.v[j];
1114
1115 build_connectivity(M_);
1116 orient_volume_mesh(M_);
1117 build_connectivity(M_);
1118
1119 M = M_;
1120 }
1121}
1122void MeshProcessing3D::refine_red_refinement_tet(Mesh3DStorage &M, int iter)
1123{
1124
1125 // double hmin=10000, hmax=0, havg=0;
1126 // for(const auto & e:M.edges){
1127 // Eigen::Vector3d v0 = M.points.col(e.vs[0]), v1 = M.points.col(e.vs[1]);
1128 // double len = (v0-v1).norm();
1129 // if(len<hmin) hmin=len;
1130 // if(len>hmax) hmax = len;
1131 // havg+=len;
1132 // }
1133 // havg/=M.edges.size();
1134 // std::cout<<"hmin, hmax, havg: "<<hmin<<" "<<hmax<<" "<<havg<<std::endl;
1135
1136 for (int i = 0; i < iter; i++)
1137 {
1138
1139 Mesh3DStorage M_;
1140 M_.type = MeshType::TET;
1141
1142 vector<int> E2V(M.edges.size());
1143
1144 int vn = 0;
1145 for (const auto &v : M.vertices)
1146 {
1147 Vertex v_;
1148 v_.id = vn++;
1149 v_.v.resize(3);
1150 for (int j = 0; j < 3; j++)
1151 v_.v[j] = M.points(j, v.id);
1152 M_.vertices.push_back(v_);
1153 }
1154
1155 for (const auto &e : M.edges)
1156 {
1157 Vertex v;
1158 v.id = vn++;
1159 v.v.resize(3);
1160
1161 Vector3d center;
1162 center.setZero();
1163 for (auto vid : e.vs)
1164 center += M.points.col(vid);
1165 center /= e.vs.size();
1166 for (int j = 0; j < 3; j++)
1167 v.v[j] = center[j];
1168
1169 M_.vertices.push_back(v);
1170 E2V[e.id] = v.id;
1171 }
1172 // for (const auto &ele : M.elements) {
1173 // Element ele_;
1174 // ele_.id = M_.elements.size();
1175 // ele_.fs.resize(4, -1);
1176 // ele_.fs_flag.resize(4, 1);
1177 // ele_.vs = ele.vs;
1178
1179 // ele_.hex = false;
1180 // ele_.v_in_Kernel.resize(3);
1181
1182 // Vector3d center;
1183 // center.setZero();
1184 // for (const auto &evid : ele_.vs) for (int j = 0; j < 3; j++)center[j] += M_.vertices[evid].v[j];
1185 // center /= ele_.vs.size();
1186 // for (int j = 0; j < 3; j++)ele_.v_in_Kernel[j] = center[j];
1187
1188 // M_.elements.push_back(ele_);
1189 // }
1190
1191 auto shared_edge = [&](int v0, int v1, int &e) -> bool {
1192 auto &es0 = M.vertices[v0].neighbor_es, &es1 = M.vertices[v1].neighbor_es;
1193 std::sort(es0.begin(), es0.end());
1194 std::sort(es1.begin(), es1.end());
1195 std::vector<uint32_t> es;
1196 std::set_intersection(es0.begin(), es0.end(), es1.begin(), es1.end(), back_inserter(es));
1197 if (es.size())
1198 {
1199 e = es[0];
1200 return true;
1201 }
1202 return false;
1203 };
1204
1205 std::vector<bool> e_flag(M.edges.size(), false);
1206
1207 for (const auto &ele : M.elements)
1208 { // 1 --> 8
1209
1210 for (short i = 0; i < 4; i++)
1211 { // four corners
1212 Element ele_;
1213 ele_.id = M_.elements.size();
1214 ele_.fs.resize(4, -1);
1215 ele_.fs_flag.resize(4, 1);
1216
1217 ele_.hex = false;
1218 ele_.v_in_Kernel.resize(3);
1219
1220 ele_.vs.push_back(ele.vs[i]);
1221 for (short j = 0; j < 4; j++)
1222 {
1223 if (j == i)
1224 continue;
1225 int v0 = ele.vs[i], v1 = ele.vs[j];
1226 int e = -1;
1227 if (shared_edge(v0, v1, e))
1228 ele_.vs.push_back(E2V[e]);
1229 }
1230 Vector3d center;
1231 center.setZero();
1232 for (const auto &evid : ele_.vs)
1233 for (int j = 0; j < 3; j++)
1234 center[j] += M_.vertices[evid].v[j];
1235 center /= ele_.vs.size();
1236 for (int j = 0; j < 3; j++)
1237 ele_.v_in_Kernel[j] = center[j];
1238
1239 M_.elements.push_back(ele_);
1240 }
1241
1242 // 6 edges
1243 std::vector<int> edges(6);
1244 for (int i = 0; i < 6; i++)
1245 {
1246 int v0 = ele.vs[tet_edges[i][0]], v1 = ele.vs[tet_edges[i][1]];
1247 int e = -1;
1248 if (shared_edge(v0, v1, e))
1249 edges[i] = e;
1250 }
1251 // the longest edge
1252 int lv0 = E2V[edges[0]], lv1 = E2V[edges[5]];
1253 e_flag[edges[0]] = true;
1254 e_flag[edges[5]] = true;
1255 for (short i = 0; i < 4; i++)
1256 { // four faces
1257 Element ele_;
1258 ele_.id = M_.elements.size();
1259 ele_.fs.resize(4, -1);
1260 ele_.fs_flag.resize(4, 1);
1261
1262 ele_.hex = false;
1263 ele_.v_in_Kernel.resize(3);
1264
1265 ele_.vs.push_back(lv0);
1266 ele_.vs.push_back(lv1);
1267 for (short j = 0; j < 3; j++)
1268 {
1269 int c_e = M.faces[ele.fs[i]].es[j];
1270 if (e_flag[c_e])
1271 continue;
1272 ele_.vs.push_back(E2V[c_e]);
1273 }
1274 Vector3d center;
1275 center.setZero();
1276 for (const auto &evid : ele_.vs)
1277 for (int j = 0; j < 3; j++)
1278 center[j] += M_.vertices[evid].v[j];
1279 center /= ele_.vs.size();
1280 for (int j = 0; j < 3; j++)
1281 ele_.v_in_Kernel[j] = center[j];
1282
1283 M_.elements.push_back(ele_);
1284 }
1285 e_flag[edges[0]] = e_flag[edges[5]] = false;
1286 }
1287
1288 M_.points.resize(3, M_.vertices.size());
1289 for (auto v : M_.vertices)
1290 for (int j = 0; j < 3; j++)
1291 M_.points(j, v.id) = v.v[j];
1292
1293 // orient tets
1294 const auto &t = M.elements[0].vs;
1295 Vector3d c0 = M.points.col(t[0]);
1296 Vector3d c1 = M.points.col(t[1]);
1297 Vector3d c2 = M.points.col(t[2]);
1298 Vector3d c3 = M.points.col(t[3]);
1299 bool signed_volume = a_jacobian(c0, c1, c2, c3) > 0 ? true : false;
1300
1301 for (auto &ele : M_.elements)
1302 {
1303 c0 = M_.points.col(ele.vs[0]);
1304 c1 = M_.points.col(ele.vs[1]);
1305 c2 = M_.points.col(ele.vs[2]);
1306 c3 = M_.points.col(ele.vs[3]);
1307 bool sign = a_jacobian(c0, c1, c2, c3) > 0 ? true : false;
1308 if (sign != signed_volume)
1309 std::swap(ele.vs[1], ele.vs[3]);
1310 }
1311
1312 // Fs
1313 std::vector<std::vector<uint32_t>> total_fs;
1314 total_fs.reserve(M.elements.size() * 4);
1315 std::vector<std::tuple<uint32_t, uint32_t, uint32_t, uint32_t, uint32_t, uint32_t>> tempF;
1316 tempF.reserve(M.elements.size() * 4);
1317 std::vector<uint32_t> vs(3);
1318 for (const auto &ele : M_.elements)
1319 {
1320 for (int i = 0; i < 4; i++)
1321 {
1322 vs[0] = ele.vs[tet_faces[i][0]];
1323 vs[1] = ele.vs[tet_faces[i][1]];
1324 vs[2] = ele.vs[tet_faces[i][2]];
1325 total_fs.push_back(vs);
1326 std::sort(vs.begin(), vs.end());
1327 tempF.push_back(std::make_tuple(vs[0], vs[1], vs[2], ele.id * 4 + i, ele.id, i));
1328 }
1329 }
1330 std::sort(tempF.begin(), tempF.end());
1331 M_.faces.reserve(tempF.size() / 3);
1332 Face f;
1333 f.boundary = true;
1334 uint32_t F_num = 0;
1335 for (uint32_t i = 0; i < tempF.size(); ++i)
1336 {
1337 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]))))
1338 {
1339 f.id = F_num;
1340 F_num++;
1341 f.vs = total_fs[std::get<3>(tempF[i])];
1342 M_.faces.push_back(f);
1343 }
1344 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])))
1345 M_.faces[F_num - 1].boundary = false;
1346
1347 M_.elements[std::get<4>(tempF[i])].fs[std::get<5>(tempF[i])] = F_num - 1;
1348 }
1349
1350 build_connectivity(M_);
1351 orient_volume_mesh(M_);
1352 build_connectivity(M_);
1353
1354 M = M_;
1355
1356 double hmin = 10000, hmax = 0, havg = 0;
1357 for (const auto &e : M.edges)
1358 {
1359 Eigen::Vector3d v0 = M.points.col(e.vs[0]), v1 = M.points.col(e.vs[1]);
1360 double len = (v0 - v1).norm();
1361 if (len < hmin)
1362 hmin = len;
1363 if (len > hmax)
1364 hmax = len;
1365 havg += len;
1366 }
1367 havg /= M.edges.size();
1368 // std::cout<<"hmin, hmax, havg: "<<hmin<<" "<<hmax<<" "<<havg<<std::endl;
1369 }
1370}
1371void MeshProcessing3D::straight_sweeping(const Mesh3DStorage &Mi, int sweep_coord, double height, int nlayer, Mesh3DStorage &Mo)
1372{
1373 if (sweep_coord > 2 || sweep_coord < 0)
1374 {
1375 logger().error("invalid sweeping direction!");
1376 return;
1377 }
1378 if (Mi.type != MeshType::H_SUR && Mi.type != MeshType::TRI && Mi.type != MeshType::QUA)
1379 {
1380 logger().error("invalid planar surface!");
1381 return;
1382 }
1383 if (height <= 0 || nlayer < 1)
1384 {
1385 logger().error("invalid height or number of layers!");
1386 return;
1387 }
1388
1389 Mo.points.resize(3, 0);
1390 Mo.vertices.clear();
1391 Mo.edges.clear();
1392 Mo.faces.clear();
1393 Mo.elements.clear();
1394 Mo.type = MeshType::HYB;
1395 // v, layers
1396 std::vector<std::vector<int>> Vlayers(nlayer + 1);
1397 std::vector<double> interval(3, 0);
1398 interval[sweep_coord] = height / nlayer;
1399
1400 for (int i = 0; i < nlayer + 1; i++)
1401 {
1402 std::vector<int> a_layer;
1403 for (auto v : Mi.vertices)
1404 {
1405 Vertex v_;
1406 v_.id = Mo.vertices.size();
1407 v_.v.resize(3);
1408 for (int j = 0; j < 3; j++)
1409 v_.v[j] = v.v[j] + i * interval[j];
1410
1411 Mo.vertices.push_back(v_);
1412 a_layer.push_back(v_.id);
1413 }
1414 Vlayers[i] = a_layer;
1415 }
1416 // f
1417 std::vector<std::vector<int>> Flayers(nlayer + 1);
1418 for (int i = 0; i < nlayer + 1; i++)
1419 {
1420 std::vector<int> a_layer;
1421 for (auto f : Mi.faces)
1422 {
1423 Face f_;
1424 f_.id = Mo.faces.size();
1425 for (auto vid : f.vs)
1426 f_.vs.push_back(Vlayers[i][vid]);
1427
1428 Mo.faces.push_back(f_);
1429 a_layer.push_back(f_.id);
1430 }
1431 Flayers[i] = a_layer;
1432 }
1433 // ef
1434 std::vector<std::vector<int>> EFlayers(nlayer);
1435 for (int i = 0; i < nlayer; i++)
1436 {
1437 std::vector<int> a_layer;
1438 for (auto e : Mi.edges)
1439 {
1440 Face f_;
1441 f_.id = Mo.faces.size();
1442 int v0 = e.vs[0], v1 = e.vs[1];
1443 f_.vs.push_back(Vlayers[i][v0]);
1444 f_.vs.push_back(Vlayers[i][v1]);
1445 f_.vs.push_back(Vlayers[i + 1][v1]);
1446 f_.vs.push_back(Vlayers[i + 1][v0]);
1447
1448 Mo.faces.push_back(f_);
1449 a_layer.push_back(f_.id);
1450 }
1451 EFlayers[i] = a_layer;
1452 }
1453 // ele
1454 for (int i = 0; i < nlayer; i++)
1455 {
1456 for (auto f : Mi.faces)
1457 {
1458 Element ele;
1459 ele.id = Mo.elements.size();
1460 ele.hex = (f.vs.size() == 4);
1461
1462 ele.fs.push_back(Flayers[i][f.id]);
1463 ele.fs.push_back(Flayers[i + 1][f.id]);
1464 for (auto eid : f.es)
1465 ele.fs.push_back(EFlayers[i][eid]);
1466
1467 ele.fs_flag.resize(ele.fs.size(), false);
1468
1469 ele.v_in_Kernel.resize(3, 0);
1470 int nv = 0;
1471 for (int j = 0; j < 2; j++)
1472 {
1473 nv += Mo.faces[ele.fs[j]].vs.size();
1474 for (auto vid : Mo.faces[ele.fs[j]].vs)
1475 for (int k = 0; k < 3; k++)
1476 ele.v_in_Kernel[k] += Mo.vertices[vid].v[k];
1477 }
1478 for (int k = 0; k < 3; k++)
1479 ele.v_in_Kernel[k] /= nv;
1480
1481 Mo.elements.push_back(ele);
1482 }
1483 }
1484
1485 Mo.points.resize(3, Mo.vertices.size());
1486 for (auto v : Mo.vertices)
1487 for (int j = 0; j < 3; j++)
1488 Mo.points(j, v.id) = v.v[j];
1489
1490 build_connectivity(Mo);
1491 orient_volume_mesh(Mo);
1492 build_connectivity(Mo);
1493}
1494
1495void MeshProcessing3D::orient_surface_mesh(Mesh3DStorage &hmi)
1496{
1497
1498 vector<bool> flag(hmi.faces.size(), true);
1499 flag[0] = false;
1500
1501 std::queue<uint32_t> pf_temp;
1502 pf_temp.push(0);
1503 while (!pf_temp.empty())
1504 {
1505 uint32_t fid = pf_temp.front();
1506 pf_temp.pop();
1507 for (auto eid : hmi.faces[fid].es)
1508 for (auto nfid : hmi.edges[eid].neighbor_fs)
1509 {
1510 if (!flag[nfid])
1511 continue;
1512 uint32_t v0 = hmi.edges[eid].vs[0], v1 = hmi.edges[eid].vs[1];
1513 int32_t v0_pos = std::find(hmi.faces[fid].vs.begin(), hmi.faces[fid].vs.end(), v0) - hmi.faces[fid].vs.begin();
1514 int32_t v1_pos = std::find(hmi.faces[fid].vs.begin(), hmi.faces[fid].vs.end(), v1) - hmi.faces[fid].vs.begin();
1515
1516 if ((v0_pos + 1) % hmi.faces[fid].vs.size() != v1_pos)
1517 swap(v0, v1);
1518
1519 int32_t v0_pos_ = std::find(hmi.faces[nfid].vs.begin(), hmi.faces[nfid].vs.end(), v0) - hmi.faces[nfid].vs.begin();
1520 int32_t v1_pos_ = std::find(hmi.faces[nfid].vs.begin(), hmi.faces[nfid].vs.end(), v1) - hmi.faces[nfid].vs.begin();
1521
1522 if ((v0_pos_ + 1) % hmi.faces[nfid].vs.size() == v1_pos_)
1523 std::reverse(hmi.faces[nfid].vs.begin(), hmi.faces[nfid].vs.end());
1524
1525 pf_temp.push(nfid);
1526 flag[nfid] = false;
1527 }
1528 }
1529 double res = 0;
1530 Vector3d ori;
1531 ori.setZero();
1532 for (auto f : hmi.faces)
1533 {
1534 auto &fvs = f.vs;
1535 Vector3d center;
1536 center.setZero();
1537 for (auto vid : fvs)
1538 center += hmi.points.col(vid);
1539 center /= fvs.size();
1540
1541 for (uint32_t j = 0; j < fvs.size(); j++)
1542 {
1543 Vector3d x = hmi.points.col(fvs[j]) - ori, y = hmi.points.col(fvs[(j + 1) % fvs.size()]) - ori, z = center - ori;
1544 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]));
1545 }
1546 }
1547 if (res > 0)
1548 {
1549 for (uint32_t i = 0; i < hmi.faces.size(); i++)
1550 std::reverse(hmi.faces[i].vs.begin(), hmi.faces[i].vs.end());
1551 }
1552}
1553void MeshProcessing3D::orient_volume_mesh(Mesh3DStorage &hmi)
1554{
1555 // surface orienting
1556 Mesh3DStorage M_sur;
1557 M_sur.type = MeshType::H_SUR;
1558 int bvn = 0;
1559 for (auto v : hmi.vertices)
1560 if (v.boundary)
1561 bvn++;
1562 M_sur.points.resize(3, bvn);
1563 bvn = 0;
1564 vector<int> V_map(hmi.vertices.size(), -1), V_map_reverse;
1565 for (auto v : hmi.vertices)
1566 if (v.boundary)
1567 {
1568 M_sur.points.col(bvn++) = hmi.points.col(v.id);
1569 Vertex v_;
1570 v_.id = M_sur.vertices.size();
1571 M_sur.vertices.push_back(v_);
1572
1573 V_map[v.id] = v_.id;
1574 V_map_reverse.push_back(v.id);
1575 }
1576 for (auto f : hmi.faces)
1577 if (f.boundary)
1578 {
1579 for (int j = 0; j < f.vs.size(); j++)
1580 {
1581 f.id = M_sur.faces.size();
1582 f.es.clear();
1583 f.neighbor_hs.clear();
1584 f.vs[j] = V_map[f.vs[j]];
1585 M_sur.vertices[f.vs[j]].neighbor_fs.push_back(f.id);
1586 }
1587 M_sur.faces.push_back(f);
1588 }
1589 build_connectivity(M_sur);
1590 orient_surface_mesh(M_sur);
1591
1592 int fn_ = 0;
1593 for (auto &f : hmi.faces)
1594 if (f.boundary)
1595 {
1596 for (int j = 0; j < f.vs.size(); j++)
1597 f.vs[j] = V_map_reverse[M_sur.faces[fn_].vs[j]];
1598 fn_++;
1599 }
1600 // volume orienting
1601 vector<bool> F_tag(hmi.faces.size(), true);
1602 std::vector<short> F_visit(hmi.faces.size(), 0); // 0 un-visited, 1 visited once, 2 visited twice
1603 for (uint32_t j = 0; j < hmi.faces.size(); j++)
1604 if (hmi.faces[j].boundary)
1605 {
1606 F_visit[j]++;
1607 }
1608 std::vector<bool> F_state(hmi.faces.size(), false); // false is the reverse direction, true is the same direction
1609 std::vector<bool> P_visit(hmi.elements.size(), false);
1610 while (true)
1611 {
1612 std::vector<uint32_t> candidates;
1613 for (uint32_t j = 0; j < F_visit.size(); j++)
1614 if (F_visit[j] == 1)
1615 candidates.push_back(j);
1616 if (!candidates.size())
1617 break;
1618 for (auto ca : candidates)
1619 {
1620 if (F_visit[ca] == 2)
1621 continue;
1622 uint32_t pid = hmi.faces[ca].neighbor_hs[0];
1623 if (P_visit[pid])
1624 if (hmi.faces[ca].neighbor_hs.size() == 2)
1625 pid = hmi.faces[ca].neighbor_hs[1];
1626 if (P_visit[pid])
1627 {
1628 logger().error("bug");
1629 }
1630 auto &fs = hmi.elements[pid].fs;
1631 for (auto fid : fs)
1632 F_tag[fid] = false;
1633
1634 uint32_t start_f = ca;
1635 F_tag[start_f] = true;
1636 F_visit[ca]++;
1637 if (F_state[ca])
1638 F_state[ca] = false;
1639 else
1640 F_state[ca] = true;
1641
1642 std::queue<uint32_t> pf_temp;
1643 pf_temp.push(start_f);
1644 while (!pf_temp.empty())
1645 {
1646 uint32_t fid = pf_temp.front();
1647 pf_temp.pop();
1648 for (auto eid : hmi.faces[fid].es)
1649 for (auto nfid : hmi.edges[eid].neighbor_fs)
1650 {
1651
1652 if (F_tag[nfid])
1653 continue;
1654 uint32_t v0 = hmi.edges[eid].vs[0], v1 = hmi.edges[eid].vs[1];
1655 int32_t v0_pos = std::find(hmi.faces[fid].vs.begin(), hmi.faces[fid].vs.end(), v0) - hmi.faces[fid].vs.begin();
1656 int32_t v1_pos = std::find(hmi.faces[fid].vs.begin(), hmi.faces[fid].vs.end(), v1) - hmi.faces[fid].vs.begin();
1657
1658 if ((v0_pos + 1) % hmi.faces[fid].vs.size() != v1_pos)
1659 std::swap(v0, v1);
1660
1661 int32_t v0_pos_ = std::find(hmi.faces[nfid].vs.begin(), hmi.faces[nfid].vs.end(), v0) - hmi.faces[nfid].vs.begin();
1662 int32_t v1_pos_ = std::find(hmi.faces[nfid].vs.begin(), hmi.faces[nfid].vs.end(), v1) - hmi.faces[nfid].vs.begin();
1663
1664 if (F_state[fid])
1665 {
1666 if ((v0_pos_ + 1) % hmi.faces[nfid].vs.size() == v1_pos_)
1667 F_state[nfid] = false;
1668 else
1669 F_state[nfid] = true;
1670 }
1671 else if (!F_state[fid])
1672 {
1673 if ((v0_pos_ + 1) % hmi.faces[nfid].vs.size() == v1_pos_)
1674 F_state[nfid] = true;
1675 else
1676 F_state[nfid] = false;
1677 }
1678
1679 F_visit[nfid]++;
1680
1681 pf_temp.push(nfid);
1682 F_tag[nfid] = true;
1683 }
1684 }
1685 P_visit[pid] = true;
1686 for (uint32_t j = 0; j < fs.size(); j++)
1687 hmi.elements[pid].fs_flag[j] = F_state[fs[j]];
1688 }
1689 }
1690}
1691void MeshProcessing3D::ele_subdivison_levels(const Mesh3DStorage &hmi, std::vector<int> &Ls)
1692{
1693
1694 Ls.clear();
1695 Ls.resize(hmi.elements.size(), 1);
1696 std::vector<double> volumes(hmi.elements.size(), 0);
1697
1698 auto compute_volume = [&](const int id, double &vol) {
1699 Vector3d ori;
1700 ori.setZero();
1701 for (auto f : hmi.elements[id].fs)
1702 {
1703 auto &fvs = hmi.faces[f].vs;
1704 Vector3d center;
1705 center.setZero();
1706 for (auto vid : fvs)
1707 center += hmi.points.col(vid);
1708 center /= fvs.size();
1709
1710 for (uint32_t j = 0; j < fvs.size(); j++)
1711 {
1712 Vector3d x = hmi.points.col(fvs[j]) - ori, y = hmi.points.col(fvs[(j + 1) % fvs.size()]) - ori, z = center - ori;
1713 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]));
1714 }
1715 }
1716 vol = std::abs(vol);
1717 };
1718
1719 for (const auto &ele : hmi.elements)
1720 compute_volume(ele.id, volumes[ele.id]);
1721
1722 double ave_volume = 0;
1723 for (const auto &v : volumes)
1724 ave_volume += v;
1725 ave_volume /= volumes.size();
1726 for (int i = 0; i < Ls.size(); i++)
1727 if (!hmi.elements[i].hex)
1728 {
1729 Ls[i] = volumes[i] / ave_volume;
1730 if (Ls[i] < 1)
1731 Ls[i] = 1;
1732 }
1733}
1734
1735// template<typename T>
1736// void MeshProcessing3D::set_intersection_own(const std::vector<T> &A, const std::vector<T> &B, std::vector<T> &C, const int &num){
1737void MeshProcessing3D::set_intersection_own(const std::vector<uint32_t> &A, const std::vector<uint32_t> &B, std::array<uint32_t, 2> &C, int &num)
1738{
1739 // void MeshProcessing3D::set_intersection_own( std::vector<uint32_t> &A, std::vector<uint32_t> &B, std::vector<uint32_t> &C, int &num)
1740 // C.resize(num);
1741 int n = 0;
1742 for (auto &a : A)
1743 {
1744 for (auto &b : B)
1745 {
1746 if (a == b)
1747 {
1748 C[n++] = a;
1749 if (n == num)
1750 break;
1751 }
1752 }
1753 if (n == num)
1754 break;
1755 }
1756}
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:71
polyfem::mesh::Mesh3DStorage::vertices
std::vector< Vertex > vertices
Definition Mesh3DStorage.hpp:68
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:679
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:1737
polyfem::mesh::MeshProcessing3D::refine_red_refinement_tet
void refine_red_refinement_tet(Mesh3DStorage &M, int iter)
Definition MeshProcessing3D.cpp:1122
polyfem::mesh::MeshProcessing3D::straight_sweeping
void straight_sweeping(const Mesh3DStorage &Mi, int sweep_coord, double height, int nlayer, Mesh3DStorage &Mo)
Definition MeshProcessing3D.cpp:1371
polyfem::mesh::MeshProcessing3D::orient_surface_mesh
void orient_surface_mesh(Mesh3DStorage &hmi)
Definition MeshProcessing3D.cpp:1495
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:1553
polyfem::mesh::MeshProcessing3D::ele_subdivison_levels
void ele_subdivison_levels(const Mesh3DStorage &hmi, std::vector< int > &Ls)
Definition MeshProcessing3D.cpp:1691
polyfem::mesh::MeshProcessing3D::scaled_jacobian
bool scaled_jacobian(Mesh3DStorage &hmi, Mesh_Quality &mq)
Definition MeshProcessing3D.cpp:608
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:714
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:681
polyfem::logger
spdlog::logger & logger()
Retrieves the current logger.
Definition Logger.cpp:42
polyfem::mesh::Edge::vs
std::vector< uint32_t > vs
Definition Mesh3DStorage.hpp:25
polyfem::mesh::Element::v_in_Kernel
std::vector< double > v_in_Kernel
Definition Mesh3DStorage.hpp:50
polyfem::mesh::Element::fs
std::vector< uint32_t > fs
Definition Mesh3DStorage.hpp:47
polyfem::mesh::Element::vs
std::vector< uint32_t > vs
Definition Mesh3DStorage.hpp:45
polyfem::mesh::Element::fs_flag
std::vector< bool > fs_flag
Definition Mesh3DStorage.hpp:48
polyfem::mesh::Face::vs
std::vector< uint32_t > vs
Definition Mesh3DStorage.hpp:35
polyfem::mesh::Mesh_Quality
Definition Mesh3DStorage.hpp:199
polyfem::mesh::Mesh_Quality::deviation_Jacobian
double deviation_Jacobian
Definition Mesh3DStorage.hpp:203
polyfem::mesh::Mesh_Quality::min_Jacobian
double min_Jacobian
Definition Mesh3DStorage.hpp:201
polyfem::mesh::Mesh_Quality::H_Js
Eigen::VectorXd H_Js
Definition Mesh3DStorage.hpp:205
polyfem::mesh::Mesh_Quality::V_Js
Eigen::VectorXd V_Js
Definition Mesh3DStorage.hpp:204
polyfem::mesh::Mesh_Quality::ave_Jacobian
double ave_Jacobian
Definition Mesh3DStorage.hpp:202
polyfem::mesh::Vertex::v
std::vector< double > v
Definition Mesh3DStorage.hpp:13