OGS
MeshRevision.cpp
Go to the documentation of this file.
1
15#include "MeshRevision.h"
16
17#include <numeric>
18#include <range/v3/algorithm/copy.hpp>
19#include <range/v3/algorithm/transform.hpp>
20#include <range/v3/range/conversion.hpp>
21#include <range/v3/view/transform.hpp>
22
23#include "BaseLib/Algorithm.h"
24#include "BaseLib/Logging.h"
25#include "GeoLib/Grid.h"
28#include "MeshLib/Mesh.h"
29#include "MeshLib/Properties.h"
31
32namespace MeshToolsLib
33{
36unsigned lutPrismThirdNode(unsigned const id1, unsigned const id2)
37{
38 if ((id1 == 0 && id2 == 1) || (id1 == 1 && id2 == 0))
39 {
40 return 2;
41 }
42 if ((id1 == 1 && id2 == 2) || (id1 == 2 && id2 == 1))
43 {
44 return 0;
45 }
46 if ((id1 == 0 && id2 == 2) || (id1 == 2 && id2 == 0))
47 {
48 return 1;
49 }
50 if ((id1 == 3 && id2 == 4) || (id1 == 4 && id2 == 3))
51 {
52 return 5;
53 }
54 if ((id1 == 4 && id2 == 5) || (id1 == 5 && id2 == 4))
55 {
56 return 3;
57 }
58 if ((id1 == 3 && id2 == 5) || (id1 == 5 && id2 == 3))
59 {
60 return 4;
61 }
62 return std::numeric_limits<unsigned>::max();
63}
64} // namespace MeshToolsLib
65
66namespace
67{
68template <typename ElementType>
69std::unique_ptr<MeshLib::Element> createElement(
70 std::span<MeshLib::Node* const> const element_nodes,
71 std::vector<MeshLib::Node*> const& nodes,
72 std::array<std::size_t, ElementType::n_all_nodes> const local_ids)
73{
74 using namespace MeshLib::views;
75 auto lookup_in = [](auto const& values)
76 {
77 return ranges::views::transform([&values](std::size_t const n)
78 { return values[n]; });
79 };
80
81 std::array<MeshLib::Node*, ElementType::n_all_nodes> new_nodes;
82 ranges::copy(local_ids | lookup_in(element_nodes) | ids | lookup_in(nodes),
83 begin(new_nodes));
84
85 return std::make_unique<ElementType>(new_nodes);
86}
87
89unsigned subdivideQuad(MeshLib::Element const* const quad,
90 std::vector<MeshLib::Node*> const& nodes,
91 std::vector<MeshLib::Element*>& new_elements)
92{
93 std::array<std::size_t, 3> const tri1_node_ids{0, 1, 2};
94 new_elements.push_back(
95 createElement<MeshLib::Tri>(quad->nodes(), nodes, tri1_node_ids)
96 .release());
97
98 std::array<std::size_t, 3> const tri2_node_ids{0, 2, 3};
99 new_elements.push_back(
100 createElement<MeshLib::Tri>(quad->nodes(), nodes, tri2_node_ids)
101 .release());
102
103 return 2;
104}
105
107unsigned subdividePrism(MeshLib::Element const* const prism,
108 std::vector<MeshLib::Node*> const& nodes,
109 std::vector<MeshLib::Element*>& new_elements)
110{
111 auto addTetrahedron =
112 [&prism, &nodes, &new_elements](std::array<std::size_t, 4> const ids)
113 {
114 new_elements.push_back(
115 createElement<MeshLib::Tet>(prism->nodes(), nodes, ids).release());
116 };
117
118 addTetrahedron({0, 1, 2, 3});
119 addTetrahedron({3, 2, 4, 5});
120 addTetrahedron({2, 1, 3, 4});
121
122 return 3;
123}
124
126unsigned subdivideHex(MeshLib::Element const* const hex,
127 std::vector<MeshLib::Node*> const& nodes,
128 std::vector<MeshLib::Element*>& new_elements)
129{
130 auto prism1 =
131 createElement<MeshLib::Prism>(hex->nodes(), nodes, {0, 2, 1, 4, 6, 5});
132 subdividePrism(prism1.get(), nodes, new_elements);
133
134 auto prism2 =
135 createElement<MeshLib::Prism>(hex->nodes(), nodes, {4, 6, 7, 0, 2, 3});
136 subdividePrism(prism2.get(), nodes, new_elements);
137
138 return 6;
139}
140
142unsigned subdividePyramid(MeshLib::Element const* const pyramid,
143 std::vector<MeshLib::Node*> const& nodes,
144 std::vector<MeshLib::Element*>& new_elements)
145{
146 auto addTetrahedron =
147 [&pyramid, &nodes, &new_elements](std::array<std::size_t, 4> const ids)
148 {
149 new_elements.push_back(
150 createElement<MeshLib::Tet>(pyramid->nodes(), nodes, ids)
151 .release());
152 };
153
154 addTetrahedron({0, 1, 2, 4});
155 addTetrahedron({0, 2, 3, 4});
156
157 return 2;
158}
159
163 const std::vector<MeshLib::Node*>& nodes)
164{
165 std::array<std::size_t, 2> line_node_ids = {0, 0};
166 for (unsigned i = 1; i < element->getNumberOfBaseNodes(); ++i)
167 {
168 if (element->getNode(i)->getID() != element->getNode(0)->getID())
169 {
170 line_node_ids[1] = i;
171 break;
172 }
173 }
174 assert(line_node_ids[1] != 0);
175 return createElement<MeshLib::Line>(element->nodes(), nodes, line_node_ids)
176 .release();
177}
178
182 const std::vector<MeshLib::Node*>& nodes)
183{
184 // TODO?
185 // In theory three unique nodes could also be reduced to two lines e.g. with
186 // a quad where two diametral nodes collapse. This case is currently not
187 // implemented!
188 std::array<MeshLib::Node*, 3> tri_nodes;
189 tri_nodes[0] = nodes[element->getNode(0)->getID()];
190 tri_nodes[2] = nullptr;
191 for (unsigned i = 1; i < element->getNumberOfBaseNodes(); ++i)
192 {
193 if (element->getNode(i)->getID() != tri_nodes[0]->getID())
194 {
195 tri_nodes[1] = nodes[element->getNode(i)->getID()];
196 for (unsigned j = i + 1; j < element->getNumberOfBaseNodes(); ++j)
197 {
198 if (element->getNode(j)->getID() != tri_nodes[1]->getID())
199 {
200 tri_nodes[2] = nodes[element->getNode(j)->getID()];
201 break;
202 }
203 }
204 if (tri_nodes[2])
205 {
206 break;
207 }
208 }
209 }
210 assert(tri_nodes[2] != nullptr);
211 return new MeshLib::Tri(tri_nodes);
212}
213
217 MeshLib::Element const* const element,
218 std::vector<MeshLib::Node*> const& nodes,
219 unsigned const min_elem_dim = 1)
220{
221 std::array<MeshLib::Node*, 4> new_nodes;
222 unsigned count(0);
223 new_nodes[count++] = nodes[element->getNode(0)->getID()];
224 for (unsigned i = 1; i < element->getNumberOfBaseNodes(); ++i)
225 {
226 if (count > 3)
227 {
228 break;
229 }
230 bool unique_node(true);
231 for (unsigned j = 0; j < i; ++j)
232 {
233 if (element->getNode(i)->getID() == element->getNode(j)->getID())
234 {
235 unique_node = false;
236 break;
237 }
238 }
239 if (unique_node)
240 {
241 new_nodes[count++] = nodes[element->getNode(i)->getID()];
242 };
243 }
244
245 // test if quad or tet
246 const bool isQuad(MathLib::isCoplanar(*new_nodes[0], *new_nodes[1],
247 *new_nodes[2], *new_nodes[3]));
248 if (isQuad && min_elem_dim < 3)
249 {
250 MeshLib::Element* elem(new MeshLib::Quad(new_nodes));
251 for (unsigned i = 1; i < 3; ++i)
252 {
253 if (elem->validate().none())
254 {
255 return elem;
256 }
257
258 // change node order if not convex
259 MeshLib::Node* tmp = new_nodes[i + 1];
260 new_nodes[i + 1] = new_nodes[i];
261 new_nodes[i] = tmp;
262 }
263 return elem;
264 }
265 if (!isQuad)
266 {
267 return new MeshLib::Tet(new_nodes);
268 }
269 // is quad but min elem dim == 3
270
271 return nullptr;
272}
273
276void reducePyramid(MeshLib::Element const* const org_elem,
277 unsigned const n_unique_nodes,
278 std::vector<MeshLib::Node*> const& nodes,
279 std::vector<MeshLib::Element*>& new_elements,
280 unsigned const min_elem_dim)
281{
282 if (n_unique_nodes == 4)
283 {
284 MeshLib::Element* elem(
285 constructFourNodeElement(org_elem, nodes, min_elem_dim));
286 if (elem)
287 {
288 new_elements.push_back(elem);
289 }
290 }
291 else if (n_unique_nodes == 3 && min_elem_dim < 3)
292 {
293 new_elements.push_back(constructTri(org_elem, nodes));
294 }
295 else if (n_unique_nodes == 2 && min_elem_dim == 1)
296 {
297 new_elements.push_back(constructLine(org_elem, nodes));
298 }
299}
300
304unsigned reducePrism(MeshLib::Element const* const org_elem,
305 unsigned const n_unique_nodes,
306 std::vector<MeshLib::Node*> const& nodes,
307 std::vector<MeshLib::Element*>& new_elements,
308 unsigned const min_elem_dim)
309{
310 auto addTetrahedron =
311 [&org_elem, &nodes, &new_elements](std::array<std::size_t, 4> const ids)
312 {
313 new_elements.push_back(
314 createElement<MeshLib::Tet>(org_elem->nodes(), nodes, ids)
315 .release());
316 };
317
318 // TODO?
319 // In theory a node from the bottom triangle and a node from the top
320 // triangle that are not connected by an edge could collapse, resulting in a
321 // combination of tri and quad elements. This case is currently not tested.
322
323 // if one of the non-triangle edges collapsed, elem can be reduced to a
324 // pyramid, otherwise it will be two tets
325 if (n_unique_nodes == 5)
326 {
327 for (unsigned i = 0; i < 5; ++i)
328 {
329 for (unsigned j = i + 1; j < 6; ++j)
330 {
331 if (i != j && org_elem->getNode(i)->getID() ==
332 org_elem->getNode(j)->getID())
333 {
334 // non triangle edge collapsed
335 if (i % 3 == j % 3)
336 {
337 addTetrahedron(
338 {(i + 1) % 3, (i + 2) % 3, i, (i + 1) % 3 + 3});
339 addTetrahedron(
340 {(i + 1) % 3 + 3, (i + 2) % 3, i, (i + 2) % 3 + 3});
341 return 2;
342 }
343
344 // triangle edge collapsed
345 const unsigned i_offset = (i > 2) ? i - 3 : i + 3;
346 const unsigned j_offset = (i > 2) ? j - 3 : j + 3;
347 const unsigned k = MeshToolsLib::lutPrismThirdNode(i, j);
348 if (k == std::numeric_limits<unsigned>::max())
349 {
350 ERR("Unexpected error during prism reduction.");
351 return 0;
352 }
353 const unsigned k_offset = (i > 2) ? k - 3 : k + 3;
354
355 addTetrahedron({i_offset, j_offset, k_offset, i});
356
357 const unsigned l =
358 (MathLib::isCoplanar(*org_elem->getNode(i_offset),
359 *org_elem->getNode(k_offset),
360 *org_elem->getNode(i),
361 *org_elem->getNode(k)))
362 ? j
363 : i;
364 const unsigned l_offset = (i > 2) ? l - 3 : l + 3;
365 addTetrahedron({l_offset, k_offset, i, k});
366 return 2;
367 }
368 }
369 }
370 }
371 else if (n_unique_nodes == 4)
372 {
373 MeshLib::Element* const elem(
374 constructFourNodeElement(org_elem, nodes, min_elem_dim));
375 if (elem)
376 {
377 new_elements.push_back(elem);
378 }
379 }
380 else if (n_unique_nodes == 3 && min_elem_dim < 3)
381 {
382 new_elements.push_back(constructTri(org_elem, nodes));
383 }
384 else if (n_unique_nodes == 2 && min_elem_dim == 1)
385 {
386 new_elements.push_back(constructLine(org_elem, nodes));
387 }
388 return 1;
389}
390
393std::array<unsigned, 4> lutHexCuttingQuadNodes(unsigned id1, unsigned id2)
394{
395 if (id1 == 0 && id2 == 1)
396 {
397 return {3, 2, 5, 4};
398 }
399 if (id1 == 1 && id2 == 2)
400 {
401 return {0, 3, 6, 5};
402 }
403 if (id1 == 2 && id2 == 3)
404 {
405 return {1, 0, 7, 6};
406 }
407 if (id1 == 3 && id2 == 0)
408 {
409 return {2, 1, 4, 7};
410 }
411 if (id1 == 4 && id2 == 5)
412 {
413 return {0, 1, 6, 7};
414 }
415 if (id1 == 5 && id2 == 6)
416 {
417 return {1, 2, 7, 4};
418 }
419 if (id1 == 6 && id2 == 7)
420 {
421 return {2, 3, 4, 5};
422 }
423 if (id1 == 7 && id2 == 4)
424 {
425 return {3, 0, 5, 6};
426 }
427 if (id1 == 0 && id2 == 4)
428 {
429 return {3, 7, 5, 1};
430 }
431 if (id1 == 1 && id2 == 5)
432 {
433 return {0, 4, 6, 2};
434 }
435 if (id1 == 2 && id2 == 6)
436 {
437 return {1, 5, 7, 3};
438 }
439 if (id1 == 3 && id2 == 7)
440 {
441 return {2, 6, 4, 0};
442 }
443 if (id1 == 1 && id2 == 0)
444 {
445 return {2, 3, 4, 5};
446 }
447 if (id1 == 2 && id2 == 1)
448 {
449 return {3, 0, 5, 6};
450 }
451 if (id1 == 3 && id2 == 2)
452 {
453 return {0, 1, 6, 7};
454 }
455 if (id1 == 0 && id2 == 3)
456 {
457 return {1, 2, 7, 4};
458 }
459 if (id1 == 5 && id2 == 4)
460 {
461 return {1, 0, 7, 6};
462 }
463 if (id1 == 6 && id2 == 5)
464 {
465 return {2, 1, 4, 7};
466 }
467 if (id1 == 7 && id2 == 6)
468 {
469 return {3, 2, 5, 4};
470 }
471 if (id1 == 4 && id2 == 7)
472 {
473 return {0, 3, 6, 5};
474 }
475 if (id1 == 4 && id2 == 0)
476 {
477 return {7, 3, 1, 5};
478 }
479 if (id1 == 5 && id2 == 1)
480 {
481 return {4, 0, 2, 6};
482 }
483 if (id1 == 6 && id2 == 2)
484 {
485 return {5, 1, 3, 7};
486 }
487 if (id1 == 7 && id2 == 3)
488 {
489 return {6, 2, 0, 4};
490 }
491
492 OGS_FATAL(
493 "lutHexCuttingQuadNodes() for nodes {} and {} does not have a valid "
494 "return value.",
495 id1, id2);
496}
497
500unsigned lutHexDiametralNode(unsigned const id)
501{
502 constexpr std::array<unsigned, 8> hex_diametral_node_ids = {
503 {6, 7, 4, 5, 2, 3, 0, 1}};
504
505 return hex_diametral_node_ids[id];
506}
507
510std::pair<unsigned, unsigned> lutHexBackNodes(unsigned const i,
511 unsigned const j,
512 unsigned const k,
513 unsigned const l)
514{
515 // collapsed edges are *not* connected
516 if (lutHexDiametralNode(i) == k)
517 {
518 return {i, lutHexDiametralNode(l)};
519 }
520 if (lutHexDiametralNode(i) == l)
521 {
522 return {i, lutHexDiametralNode(k)};
523 }
524 if (lutHexDiametralNode(j) == k)
525 {
526 return {j, lutHexDiametralNode(l)};
527 }
528 if (lutHexDiametralNode(j) == l)
529 {
530 return {j, lutHexDiametralNode(k)};
531 }
532
533 // collapsed edges *are* connected
534 if (i == k)
535 {
536 return {lutHexDiametralNode(l), j};
537 }
538 if (i == l)
539 {
540 return {lutHexDiametralNode(k), j};
541 }
542 if (j == k)
543 {
544 return {lutHexDiametralNode(l), i};
545 }
546 if (j == l)
547 {
548 return {lutHexDiametralNode(k), i};
549 }
550
551 return {std::numeric_limits<unsigned>::max(),
552 std::numeric_limits<unsigned>::max()};
553}
554
555// In an element with 5 unique nodes, return the node that will be the top of
556// the resulting pyramid.
557unsigned findPyramidTopNode(MeshLib::Element const& element,
558 std::array<std::size_t, 4> const& base_node_ids)
559{
560 const std::size_t nNodes(element.getNumberOfBaseNodes());
561 for (std::size_t i = 0; i < nNodes; ++i)
562 {
563 bool top_node = true;
564 for (unsigned j = 0; j < 4; ++j)
565 {
566 if (element.getNode(i)->getID() == base_node_ids[j])
567 {
568 top_node = false;
569 }
570 }
571 if (top_node)
572 {
573 return i;
574 }
575 }
576 return std::numeric_limits<unsigned>::max(); // should never be reached if
577 // called correctly
578}
579
583unsigned reduceHex(MeshLib::Element const* const org_elem,
584 unsigned const n_unique_nodes,
585 std::vector<MeshLib::Node*> const& nodes,
586 std::vector<MeshLib::Element*>& new_elements,
587 unsigned const min_elem_dim)
588{
589 // TODO?
590 // if two diametral nodes collapse, all kinds of bizarre (2D-)element
591 // combinations could be the result. this case is currently not implemented!
592
593 if (n_unique_nodes == 7)
594 {
595 // reduce to prism + pyramid
596 for (unsigned i = 0; i < 7; ++i)
597 {
598 for (unsigned j = i + 1; j < 8; ++j)
599 {
600 if (org_elem->getNode(i)->getID() ==
601 org_elem->getNode(j)->getID())
602 {
603 const std::array<unsigned, 4> base_node_ids(
605 std::array<std::size_t, 5> const pyr_node_ids = {
606 base_node_ids[0], base_node_ids[1], base_node_ids[2],
607 base_node_ids[3], i};
608 new_elements.push_back(
610 nodes, pyr_node_ids)
611 .release());
612
613 if (i < 4 && j >= 4)
614 {
615 std::swap(i, j);
616 }
617 std::array<std::size_t, 6> const prism_node_ids{
618 base_node_ids[0], base_node_ids[3],
619 lutHexDiametralNode(j), base_node_ids[1],
620 base_node_ids[2], lutHexDiametralNode(i)};
621 new_elements.push_back(
622 createElement<MeshLib::Prism>(org_elem->nodes(), nodes,
623 prism_node_ids)
624 .release());
625 return 2;
626 }
627 }
628 }
629 }
630 else if (n_unique_nodes == 6)
631 {
632 // reduce to prism
633 for (unsigned i = 0; i < 6; ++i)
634 {
635 const MeshLib::Element* face(org_elem->getFace(i));
636 if (face->getNode(0)->getID() == face->getNode(1)->getID() &&
637 face->getNode(2)->getID() == face->getNode(3)->getID())
638 {
639 std::array<std::size_t, 6> const prism_node_ids{
641 getNodeIDinElement(*org_elem, face->getNode(0))),
643 getNodeIDinElement(*org_elem, face->getNode(1))),
644 getNodeIDinElement(*org_elem, face->getNode(2)),
646 getNodeIDinElement(*org_elem, face->getNode(2))),
648 getNodeIDinElement(*org_elem, face->getNode(3))),
649 getNodeIDinElement(*org_elem, face->getNode(0))};
650
651 new_elements.push_back(
652 createElement<MeshLib::Prism>(org_elem->nodes(), nodes,
653 prism_node_ids)
654 .release());
655 delete face;
656 return 1;
657 }
658 if (face->getNode(0)->getID() == face->getNode(3)->getID() &&
659 face->getNode(1)->getID() == face->getNode(2)->getID())
660 {
661 std::array<std::size_t, 6> const prism_node_ids{
663 getNodeIDinElement(*org_elem, face->getNode(0))),
665 getNodeIDinElement(*org_elem, face->getNode(3))),
666 getNodeIDinElement(*org_elem, face->getNode(2)),
668 getNodeIDinElement(*org_elem, face->getNode(1))),
670 getNodeIDinElement(*org_elem, face->getNode(2))),
671 getNodeIDinElement(*org_elem, face->getNode(0))};
672 new_elements.push_back(
673 createElement<MeshLib::Prism>(org_elem->nodes(), nodes,
674 prism_node_ids)
675 .release());
676 delete face;
677 return 1;
678 }
679 delete face;
680 }
681 // reduce to four tets -> divide into 2 prisms such that each has one
682 // collapsed node
683 for (unsigned i = 0; i < 7; ++i)
684 {
685 for (unsigned j = i + 1; j < 8; ++j)
686 {
687 if (org_elem->getNode(i)->getID() ==
688 org_elem->getNode(j)->getID())
689 {
690 for (unsigned k = i; k < 7; ++k)
691 {
692 for (unsigned l = k + 1; l < 8; ++l)
693 {
694 if (!(i == k && j == l) && org_elem->isEdge(i, j) &&
695 org_elem->isEdge(k, l) &&
696 org_elem->getNode(k)->getID() ==
697 org_elem->getNode(l)->getID())
698 {
699 const std::pair<unsigned, unsigned> back(
700 lutHexBackNodes(i, j, k, l));
701 if (back.first ==
702 std::numeric_limits<unsigned>::max() ||
703 back.second ==
704 std::numeric_limits<unsigned>::max())
705 {
706 ERR("Unexpected error during Hex "
707 "reduction");
708 return 0;
709 }
710
711 std::array<unsigned, 4> const cutting_plane(
712 lutHexCuttingQuadNodes(back.first,
713 back.second));
714 std::array<std::size_t, 6> const pris1_node_ids{
715 back.first, cutting_plane[0],
716 cutting_plane[3], back.second,
717 cutting_plane[1], cutting_plane[2]};
718 auto prism1 = createElement<MeshLib::Prism>(
719 org_elem->nodes(), nodes, pris1_node_ids);
720 unsigned nNewElements =
721 reducePrism(prism1.get(), 5, nodes,
722 new_elements, min_elem_dim);
723
724 std::array<std::size_t, 6> const pris2_node_ids{
725 lutHexDiametralNode(back.first),
726 cutting_plane[0],
727 cutting_plane[3],
728 lutHexDiametralNode(back.second),
729 cutting_plane[1],
730 cutting_plane[2]};
731 auto prism2 = createElement<MeshLib::Prism>(
732 org_elem->nodes(), nodes, pris2_node_ids);
733 nNewElements +=
734 reducePrism(prism2.get(), 5, nodes,
735 new_elements, min_elem_dim);
736 return nNewElements;
737 }
738 }
739 }
740 }
741 }
742 }
743 }
744 else if (n_unique_nodes == 5)
745 {
746 MeshLib::Element* tet1(constructFourNodeElement(org_elem, nodes));
747 std::array<std::size_t, 4> const first_four_node_ids = {
748 {tet1->getNode(0)->getID(), tet1->getNode(1)->getID(),
749 tet1->getNode(2)->getID(), tet1->getNode(3)->getID()}};
750 unsigned const fifth_node =
751 findPyramidTopNode(*org_elem, first_four_node_ids);
752
753 bool tet_changed(false);
755 {
756 delete tet1;
757 tet_changed = true;
758 std::array const tet1_nodes = {
759 nodes[first_four_node_ids[0]], nodes[first_four_node_ids[1]],
760 nodes[first_four_node_ids[2]],
761 nodes[org_elem->getNode(fifth_node)->getID()]};
762 new_elements.push_back(new MeshLib::Tet(tet1_nodes));
763 }
764 else
765 {
766 new_elements.push_back(tet1);
767 }
768
769 std::array const tet2_nodes = {
770 (tet_changed) ? nodes[first_four_node_ids[0]]
771 : nodes[first_four_node_ids[1]],
772 nodes[first_four_node_ids[2]], nodes[first_four_node_ids[3]],
773 nodes[org_elem->getNode(fifth_node)->getID()]};
774 new_elements.push_back(new MeshLib::Tet(tet2_nodes));
775 return 2;
776 }
777 else if (n_unique_nodes == 4)
778 {
779 MeshLib::Element* elem(
780 constructFourNodeElement(org_elem, nodes, min_elem_dim));
781 if (elem)
782 {
783 new_elements.push_back(elem);
784 return 1;
785 }
786 }
787 else if (n_unique_nodes == 3 && min_elem_dim < 3)
788 {
789 new_elements.push_back(constructTri(org_elem, nodes));
790 return 1;
791 }
792 else if (min_elem_dim == 1)
793 {
794 new_elements.push_back(constructLine(org_elem, nodes));
795 return 1;
796 }
797 return 0;
798}
799
807std::size_t subdivideElement(MeshLib::Element const* const element,
808 std::vector<MeshLib::Node*> const& nodes,
809 std::vector<MeshLib::Element*>& elements)
810{
812 {
813 return subdivideQuad(element, nodes, elements);
814 }
816 {
817 return subdivideHex(element, nodes, elements);
818 }
820 {
821 return subdividePyramid(element, nodes, elements);
822 }
824 {
825 return subdividePrism(element, nodes, elements);
826 }
827 return 0;
828}
829
830// Revises an element by removing collapsed nodes, using the nodes vector from
831// the result mesh.
832std::size_t reduceElement(MeshLib::Element const* const element,
833 unsigned const n_unique_nodes,
834 std::vector<MeshLib::Node*> const& nodes,
835 std::vector<MeshLib::Element*>& elements,
836 unsigned const min_elem_dim)
837{
838 /***************
839 * TODO: modify neighbouring elements if one elements has been subdivided
840 ***************/
842 min_elem_dim == 1)
843 {
844 elements.push_back(constructLine(element, nodes));
845 return 1;
846 }
847 if ((element->getGeomType() == MeshLib::MeshElemType::QUAD) ||
849 {
850 if (n_unique_nodes == 3 && min_elem_dim < 3)
851 {
852 elements.push_back(constructTri(element, nodes));
853 }
854 else if (min_elem_dim == 1)
855 {
856 elements.push_back(constructLine(element, nodes));
857 }
858 return 1;
859 }
861 {
862 return reduceHex(element, n_unique_nodes, nodes, elements,
863 min_elem_dim);
864 }
866 {
867 reducePyramid(element, n_unique_nodes, nodes, elements, min_elem_dim);
868 return 1;
869 }
871 {
872 return reducePrism(element, n_unique_nodes, nodes, elements,
873 min_elem_dim);
874 }
875
876 ERR("Unknown element type.");
877 return 0;
878}
881unsigned getNumberOfUniqueNodes(MeshLib::Element const* const element)
882{
883 unsigned const nNodes(element->getNumberOfBaseNodes());
884 unsigned count(nNodes);
885
886 for (unsigned i = 0; i < nNodes - 1; ++i)
887 {
888 for (unsigned j = i + 1; j < nNodes; ++j)
889 {
890 if (element->getNode(i)->getID() == element->getNode(j)->getID())
891 {
892 count--;
893 break;
894 }
895 }
896 }
897 return count;
898}
899
900template <typename T>
902 MeshLib::PropertyVector<T> const& old_prop,
903 std::vector<size_t> const& node_ids)
904{
905 std::size_t const n_nodes = node_ids.size();
906 for (std::size_t i = 0; i < n_nodes; ++i)
907 {
908 if (node_ids[i] != i)
909 {
910 continue;
911 }
912 new_prop.push_back(old_prop[i]);
913 }
914}
915
916template <typename T>
918 MeshLib::PropertyVector<T> const& old_prop,
919 std::vector<size_t> const& elem_ids)
920{
921 assert(new_prop.size() == old_prop.size());
922 assert(new_prop.size() == elem_ids.size());
923 ranges::transform(elem_ids, new_prop.begin(),
924 [&](std::size_t const i) { return old_prop[i]; });
925}
926
931 std::vector<std::size_t> const& node_ids,
932 std::vector<std::size_t> const& elem_ids)
933{
934 auto const prop_names = props.getPropertyVectorNames();
935 MeshLib::Properties new_properties;
936
937 for (auto name : prop_names)
938 {
940 1))
941 {
942 auto const* p = props.getPropertyVector<int>(
944 auto new_node_vec = new_properties.createNewPropertyVector<int>(
946 fillNodeProperty(*new_node_vec, *p, node_ids);
947 continue;
948 }
950 1))
951 {
952 auto const* p = props.getPropertyVector<float>(
954 auto new_node_vec = new_properties.createNewPropertyVector<float>(
956 fillNodeProperty(*new_node_vec, *p, node_ids);
957 continue;
958 }
959 if (props.existsPropertyVector<double>(name,
961 {
962 auto const* p = props.getPropertyVector<double>(
964 auto new_node_vec = new_properties.createNewPropertyVector<double>(
966 fillNodeProperty(*new_node_vec, *p, node_ids);
967 continue;
968 }
970 1))
971 {
972 auto const* p = props.getPropertyVector<int>(
974 auto new_cell_vec = new_properties.createNewPropertyVector<int>(
975 name, MeshLib::MeshItemType::Cell, p->size(), 1);
976 fillElemProperty(*new_cell_vec, *p, elem_ids);
977 continue;
978 }
980 1))
981 {
982 auto const* p = props.getPropertyVector<float>(
984 auto new_cell_vec = new_properties.createNewPropertyVector<float>(
985 name, MeshLib::MeshItemType::Cell, p->size(), 1);
986 fillElemProperty(*new_cell_vec, *p, elem_ids);
987 continue;
988 }
989 if (props.existsPropertyVector<double>(name,
991 {
992 auto const* p = props.getPropertyVector<double>(
994 auto new_cell_vec = new_properties.createNewPropertyVector<double>(
995 name, MeshLib::MeshItemType::Cell, p->size(), 1);
996 fillElemProperty(*new_cell_vec, *p, elem_ids);
997 continue;
998 }
999 WARN("PropertyVector {:s} not being converted.", name);
1000 }
1001 return new_properties;
1002}
1003} // namespace
1004
1005namespace MeshToolsLib
1006{
1008
1009unsigned MeshRevision::getNumberOfCollapsibleNodes(double const eps) const
1010{
1011 std::vector<std::size_t> const id_map = collapseNodeIndices(eps);
1012 std::size_t const nNodes = id_map.size();
1013 unsigned count(0);
1014 for (std::size_t i = 0; i < nNodes; ++i)
1015 {
1016 if (i != id_map[i])
1017 {
1018 count++;
1019 }
1020 }
1021 return count;
1022}
1023
1024MeshLib::Mesh* MeshRevision::simplifyMesh(const std::string& new_mesh_name,
1025 double const eps,
1026 unsigned const min_elem_dim) const
1027{
1028 if (this->_mesh.getNumberOfElements() == 0)
1029 {
1030 return nullptr;
1031 }
1032
1033 std::vector<MeshLib::Element*> const& elements(this->_mesh.getElements());
1034 auto const node_ids = collapseNodeIndices(eps);
1035 std::vector<MeshLib::Node*> new_nodes =
1036 this->constructNewNodesArray(node_ids);
1037 std::vector<MeshLib::Element*> new_elements;
1038 std::vector<std::size_t> element_ids;
1039
1040 for (std::size_t k(0); k < elements.size(); ++k)
1041 {
1042 MeshLib::Element const* const elem(elements[k]);
1043 unsigned const n_unique_nodes(getNumberOfUniqueNodes(elem));
1044 if (n_unique_nodes == elem->getNumberOfBaseNodes() &&
1045 elem->getDimension() >= min_elem_dim)
1046 {
1047 ElementErrorCode const e = elem->validate();
1049 {
1050 std::size_t const n_new_elements(
1051 subdivideElement(elem, new_nodes, new_elements));
1052 if (n_new_elements == 0)
1053 {
1054 ERR("Element {:d} has unknown element type.", k);
1056 BaseLib::cleanupVectorElements(new_nodes, new_elements);
1057 return nullptr;
1058 }
1059 element_ids.insert(element_ids.end(), n_new_elements, k);
1060 }
1061 else
1062 {
1063 new_elements.push_back(MeshLib::copyElement(elem, new_nodes));
1064 element_ids.push_back(k);
1065 }
1066 }
1067 else if (n_unique_nodes < elem->getNumberOfBaseNodes() &&
1068 n_unique_nodes > 1)
1069 {
1070 std::size_t const n_new_elements(reduceElement(
1071 elem, n_unique_nodes, new_nodes, new_elements, min_elem_dim));
1072 element_ids.insert(element_ids.end(), n_new_elements, k);
1073 }
1074 else
1075 {
1076 ERR("Something is wrong, more unique nodes than actual nodes");
1077 }
1078 }
1079
1080 auto const& props = _mesh.getProperties();
1081 MeshLib::Properties const new_properties =
1082 copyProperties(props, node_ids, element_ids);
1083
1085 if (!new_elements.empty())
1086 {
1087 return new MeshLib::Mesh(new_mesh_name, new_nodes, new_elements,
1088 true /* compute_element_neighbors */,
1089 new_properties);
1090 }
1091
1092 BaseLib::cleanupVectorElements(new_nodes, new_elements);
1093 return nullptr;
1094}
1095
1096std::vector<std::size_t> MeshRevision::collapseNodeIndices(
1097 double const eps) const
1098{
1099 const std::vector<MeshLib::Node*>& nodes(_mesh.getNodes());
1100 const std::size_t nNodes(_mesh.getNumberOfNodes());
1101 std::vector<std::size_t> id_map(nNodes);
1102 const double half_eps(eps / 2.0);
1103 const double sqr_eps(eps * eps);
1104 std::iota(id_map.begin(), id_map.end(), 0);
1105
1106 GeoLib::Grid<MeshLib::Node> const grid(nodes.begin(), nodes.end(), 64);
1107
1108 for (std::size_t k = 0; k < nNodes; ++k)
1109 {
1110 MeshLib::Node const* const node(nodes[k]);
1111 if (node->getID() != k)
1112 {
1113 continue;
1114 }
1115 std::vector<std::vector<MeshLib::Node*> const*> const node_vectors(
1116 grid.getPntVecsOfGridCellsIntersectingCube(*node, half_eps));
1117
1118 const std::size_t nVectors(node_vectors.size());
1119 for (std::size_t i = 0; i < nVectors; ++i)
1120 {
1121 const std::vector<MeshLib::Node*>& cell_vector(*node_vectors[i]);
1122 const std::size_t nGridCellNodes(cell_vector.size());
1123 for (std::size_t j = 0; j < nGridCellNodes; ++j)
1124 {
1125 MeshLib::Node const* const test_node(cell_vector[j]);
1126 // are node indices already identical (i.e. nodes will be
1127 // collapsed)
1128 if (id_map[node->getID()] == id_map[test_node->getID()])
1129 {
1130 continue;
1131 }
1132
1133 // if test_node has already been collapsed to another node x,
1134 // ignore it (if the current node would need to be collapsed
1135 // with x it would already have happened when x was tested)
1136 if (test_node->getID() != id_map[test_node->getID()])
1137 {
1138 continue;
1139 }
1140
1141 // calc distance
1142 if (MathLib::sqrDist(*node, *test_node) < sqr_eps)
1143 {
1144 id_map[test_node->getID()] = node->getID();
1145 }
1146 }
1147 }
1148 }
1149 return id_map;
1150}
1151
1152std::vector<MeshLib::Node*> MeshRevision::constructNewNodesArray(
1153 const std::vector<std::size_t>& id_map) const
1154{
1155 const std::vector<MeshLib::Node*>& nodes(_mesh.getNodes());
1156 const std::size_t nNodes(nodes.size());
1157 std::vector<MeshLib::Node*> new_nodes;
1158 new_nodes.reserve(nNodes);
1159 for (std::size_t k = 0; k < nNodes; ++k)
1160 {
1161 // all nodes that have not been collapsed with other nodes are copied
1162 // into new array
1163 if (nodes[k]->getID() == id_map[k])
1164 {
1165 std::size_t const id(new_nodes.size());
1166 new_nodes.push_back(new MeshLib::Node(
1167 (*nodes[k])[0], (*nodes[k])[1], (*nodes[k])[2], id));
1168 nodes[k]->setID(id); // the node in the old array gets the index of
1169 // the same node in the new array
1170 }
1171 // the other nodes are not copied and get the index of the nodes they
1172 // will have been collapsed with
1173 else
1174 {
1175 nodes[k]->setID(nodes[id_map[k]]->getID());
1176 }
1177 }
1178 return new_nodes;
1179}
1180
1181} // namespace MeshToolsLib
Definition of Duplicate functions.
#define OGS_FATAL(...)
Definition Error.h:26
Definition of the Grid class.
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:45
void WARN(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:40
Definition of the class Properties that implements a container of properties.
Definition of the MeshRevision class.
Definition of the Mesh class.
Collects error flags for mesh elements.
std::vector< std::vector< POINT * > const * > getPntVecsOfGridCellsIntersectingCube(P const &center, double half_len) const
Definition Grid.h:254
std::size_t getID() const
virtual MeshElemType getGeomType() const =0
virtual ElementErrorCode validate() const =0
virtual const Element * getFace(unsigned i) const =0
Returns the i-th face of the element.
virtual unsigned getNumberOfBaseNodes() const =0
virtual const Node * getNode(unsigned idx) const =0
virtual bool isEdge(unsigned i, unsigned j) const =0
Returns true if these two indices form an edge and false otherwise.
virtual constexpr unsigned getDimension() const =0
Get dimension of the mesh element.
constexpr std::span< Node *const > nodes() const
Span of element's nodes, their pointers actually.
Definition Element.h:74
std::vector< Node * > const & getNodes() const
Get the nodes-vector for the mesh.
Definition Mesh.h:108
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition Mesh.h:111
void resetNodeIDs()
Resets the IDs of all mesh-nodes to their position in the node vector.
Definition Mesh.cpp:160
Properties & getProperties()
Definition Mesh.h:136
std::size_t getNumberOfNodes() const
Get the number of nodes.
Definition Mesh.h:102
std::size_t getNumberOfElements() const
Get the number of elements.
Definition Mesh.h:99
Property manager on mesh items. Class Properties manages scalar, vector or matrix properties....
Definition Properties.h:33
std::vector< std::string > getPropertyVectorNames() const
bool existsPropertyVector(std::string_view name) const
Definition Properties.h:93
PropertyVector< T > * createNewPropertyVector(std::string_view name, MeshItemType mesh_item_type, std::size_t n_components=1)
Definition Properties.h:17
PropertyVector< T > const * getPropertyVector(std::string_view name) const
constexpr std::size_t size() const
constexpr void push_back(const PROP_VAL_TYPE &value)
constexpr PROP_VAL_TYPE * begin()
MeshRevision(MeshLib::Mesh &mesh)
unsigned getNumberOfCollapsibleNodes(double eps=std::numeric_limits< double >::epsilon()) const
Returns the number of potentially collapsible nodes.
MeshLib::Mesh & _mesh
The original mesh used for constructing the class.
MeshLib::Mesh * simplifyMesh(const std::string &new_mesh_name, double eps, unsigned min_elem_dim=1) const
std::vector< MeshLib::Node * > constructNewNodesArray(const std::vector< std::size_t > &id_map) const
std::vector< std::size_t > collapseNodeIndices(double eps) const
void cleanupVectorElements(std::vector< T * > &items)
Definition Algorithm.h:256
bool isCoplanar(const MathLib::Point3d &a, const MathLib::Point3d &b, const MathLib::Point3d &c, const MathLib::Point3d &d)
Checks if the four given points are located on a plane.
double sqrDist(MathLib::Point3d const &p0, MathLib::Point3d const &p1)
Definition Point3d.cpp:26
MeshLib specific, lazy, non-owning, non-mutating, composable range views.
Definition Mesh.h:225
TemplateElement< MeshLib::TetRule4 > Tet
Definition Tet.h:25
TemplateElement< MeshLib::TriRule3 > Tri
Definition Tri.h:26
Element * copyElement(Element const *const element, const std::vector< Node * > &nodes, std::vector< std::size_t > const *const id_map)
unsigned lutPrismThirdNode(unsigned const id1, unsigned const id2)
std::size_t reduceElement(MeshLib::Element const *const element, unsigned const n_unique_nodes, std::vector< MeshLib::Node * > const &nodes, std::vector< MeshLib::Element * > &elements, unsigned const min_elem_dim)
unsigned lutHexDiametralNode(unsigned const id)
MeshLib::Element * constructTri(MeshLib::Element const *const element, const std::vector< MeshLib::Node * > &nodes)
MeshLib::Element * constructLine(MeshLib::Element const *const element, const std::vector< MeshLib::Node * > &nodes)
unsigned reduceHex(MeshLib::Element const *const org_elem, unsigned const n_unique_nodes, std::vector< MeshLib::Node * > const &nodes, std::vector< MeshLib::Element * > &new_elements, unsigned const min_elem_dim)
unsigned reducePrism(MeshLib::Element const *const org_elem, unsigned const n_unique_nodes, std::vector< MeshLib::Node * > const &nodes, std::vector< MeshLib::Element * > &new_elements, unsigned const min_elem_dim)
void fillNodeProperty(MeshLib::PropertyVector< T > &new_prop, MeshLib::PropertyVector< T > const &old_prop, std::vector< size_t > const &node_ids)
unsigned getNumberOfUniqueNodes(MeshLib::Element const *const element)
unsigned subdivideHex(MeshLib::Element const *const hex, std::vector< MeshLib::Node * > const &nodes, std::vector< MeshLib::Element * > &new_elements)
Subdivides a Hex with nonplanar faces into tets.
std::unique_ptr< MeshLib::Element > createElement(std::span< MeshLib::Node *const > const element_nodes, std::vector< MeshLib::Node * > const &nodes, std::array< std::size_t, ElementType::n_all_nodes > const local_ids)
std::size_t subdivideElement(MeshLib::Element const *const element, std::vector< MeshLib::Node * > const &nodes, std::vector< MeshLib::Element * > &elements)
void reducePyramid(MeshLib::Element const *const org_elem, unsigned const n_unique_nodes, std::vector< MeshLib::Node * > const &nodes, std::vector< MeshLib::Element * > &new_elements, unsigned const min_elem_dim)
MeshLib::Properties copyProperties(MeshLib::Properties const &props, std::vector< std::size_t > const &node_ids, std::vector< std::size_t > const &elem_ids)
unsigned subdividePrism(MeshLib::Element const *const prism, std::vector< MeshLib::Node * > const &nodes, std::vector< MeshLib::Element * > &new_elements)
Subdivides a prism with nonplanar quad faces into two tets.
std::array< unsigned, 4 > lutHexCuttingQuadNodes(unsigned id1, unsigned id2)
std::pair< unsigned, unsigned > lutHexBackNodes(unsigned const i, unsigned const j, unsigned const k, unsigned const l)
unsigned subdividePyramid(MeshLib::Element const *const pyramid, std::vector< MeshLib::Node * > const &nodes, std::vector< MeshLib::Element * > &new_elements)
Subdivides a pyramid with a nonplanar base into two tets.
void fillElemProperty(MeshLib::PropertyVector< T > &new_prop, MeshLib::PropertyVector< T > const &old_prop, std::vector< size_t > const &elem_ids)
unsigned findPyramidTopNode(MeshLib::Element const &element, std::array< std::size_t, 4 > const &base_node_ids)
MeshLib::Element * constructFourNodeElement(MeshLib::Element const *const element, std::vector< MeshLib::Node * > const &nodes, unsigned const min_elem_dim=1)
unsigned subdivideQuad(MeshLib::Element const *const quad, std::vector< MeshLib::Node * > const &nodes, std::vector< MeshLib::Element * > &new_elements)
Subdivides a nonplanar quad into two triangles.