deal.II version GIT relicensing-2645-g9d42f38b55 2025-02-15 04:20:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
connectivity.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2020 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_tria_connectivity_h
16#define dealii_tria_connectivity_h
17
18#include <deal.II/base/config.h>
19
22
26
27#include <numeric>
28
29
31
32
33namespace internal
34{
35 namespace TriangulationImplementation
36 {
42 {
46 virtual ~CellTypeBase() = default;
47
51 virtual unsigned int
52 n_entities(const unsigned int d) const
53 {
55 (void)d;
56
57 return 0;
58 }
59
63 virtual ::ArrayView<const unsigned int>
64 vertices_of_entity(const unsigned int d, const unsigned int e) const
65 {
67 (void)d;
68 (void)e;
69
70 return {};
71 }
72
76 virtual ReferenceCell
77 type_of_entity(const unsigned int d, const unsigned int e) const
78 {
80 (void)d;
81 (void)e;
82
84 }
85
89 virtual unsigned int
90 n_lines_of_surface(const unsigned int face) const
91 {
93 (void)face;
94
95 return 0;
96 }
97
101 virtual unsigned int
102 nth_line_of_surface(const unsigned int line,
103 const unsigned int face) const
104 {
106 (void)line;
107 (void)face;
108
109 return 0;
110 }
111
115 virtual const std::array<unsigned int, 2> &
116 vertices_of_nth_line_of_surface(const unsigned int line,
117 const unsigned int face) const
118 {
120 (void)line;
121 (void)face;
122
123 const static std::array<unsigned int, 2> table = {};
124
125 return table;
126 }
127 };
128
129
130
135 {
137 vertices_of_entity(const unsigned int d,
138 const unsigned int e) const override
139 {
140 (void)e;
141
142 if (d == 1)
143 {
144 static const std::array<unsigned int, 2> table = {{0, 1}};
145
146 AssertDimension(e, 0);
147
148 return {table};
149 }
150
152
153 return {};
154 }
155
157 type_of_entity(const unsigned int d, const unsigned int e) const override
158 {
159 (void)e;
160
161 if (d == 1)
163
165
167 }
168
169 unsigned int
170 n_entities(const unsigned int d) const override
171 {
172 static std::array<unsigned int, 2> table = {{2, 1}};
173 return table[d];
174 }
175 };
176
177
178
183 {
185 vertices_of_entity(const unsigned int d,
186 const unsigned int e) const override
187 {
188 if (d == 2)
189 {
190 static const std::array<unsigned int, 3> table = {{0, 1, 2}};
191
192 AssertDimension(e, 0);
193
194 return {table};
195 }
196
197 if (d == 1)
198 {
199 static const ::ndarray<unsigned int, 3, 2> table = {
200 {{{0, 1}}, {{1, 2}}, {{2, 0}}}};
201
202 return {table[e]};
203 }
204
206
207 return {};
208 }
209
210 virtual ReferenceCell
211 type_of_entity(const unsigned int d, const unsigned int e) const override
212 {
213 (void)e;
214
215 if (d == 2)
217
218 if (d == 1)
220
222
224 }
225
226 unsigned int
227 n_entities(const unsigned int d) const override
228 {
229 static std::array<unsigned int, 3> table = {{3, 3, 1}};
230 return table[d];
231 }
232 };
233
234
235
240 {
242 vertices_of_entity(const unsigned int d,
243 const unsigned int e) const override
244 {
245 if (d == 2)
246 {
247 static const std::array<unsigned int, 4> table = {{0, 1, 2, 3}};
248
249 AssertDimension(e, 0);
250
251 return {table};
252 }
253
254 if (d == 1)
255 {
256 static const ::ndarray<unsigned int, 4, 2> table = {
257 {{{0, 2}}, {{1, 3}}, {{0, 1}}, {{2, 3}}}};
258
259 return {table[e]};
260 }
261
263
264 return {};
265 }
266
267 virtual ReferenceCell
268 type_of_entity(const unsigned int d, const unsigned int e) const override
269 {
270 (void)e;
271
272 if (d == 2)
274
275 if (d == 1)
277
279
281 }
282
283 unsigned int
284 n_entities(const unsigned int d) const override
285 {
286 static std::array<unsigned int, 3> table = {{4, 4, 1}};
287 return table[d];
288 }
289 };
290
291
292
297 {
299 vertices_of_entity(const unsigned int d,
300 const unsigned int e) const override
301 {
302 if (d == 3)
303 {
304 static const std::array<unsigned int, 4> table = {{0, 1, 2, 3}};
305
306 AssertDimension(e, 0);
307
308 return {table};
309 }
310
311 if (d == 2)
312 {
313 static const ::ndarray<unsigned int, 4, 3> table = {
314 {{{0, 1, 2}}, {{1, 0, 3}}, {{0, 2, 3}}, {{2, 1, 3}}}};
315
316 return {table[e]};
317 }
318
319 if (d == 1)
320 {
321 static const ::ndarray<unsigned int, 6, 2> table = {
322 {{{0, 1}}, {{1, 2}}, {{2, 0}}, {{0, 3}}, {{1, 3}}, {{2, 3}}}};
323
324 return {table[e]};
325 }
326
328
329 return {};
330 }
331
332 virtual ReferenceCell
333 type_of_entity(const unsigned int d, const unsigned int e) const override
334 {
335 (void)e;
336
337 if (d == 3)
339
340 if (d == 2)
342
343 if (d == 1)
345
347
349 }
350
351 unsigned int
352 n_entities(const unsigned int d) const override
353 {
354 static std::array<unsigned int, 4> table = {{4, 6, 4, 1}};
355 return table[d];
356 }
357
358 unsigned int
359 n_lines_of_surface(const unsigned int line) const override
360 {
361 (void)line;
362 return 3;
363 }
364
365 unsigned int
366 nth_line_of_surface(const unsigned int line,
367 const unsigned int face) const override
368 {
369 const static ::ndarray<unsigned int, 4, 3> table = {
370 {{{0, 1, 2}}, {{0, 3, 4}}, {{2, 5, 3}}, {{1, 4, 5}}}};
371
372 return table[face][line];
373 }
374
375 const std::array<unsigned int, 2> &
376 vertices_of_nth_line_of_surface(const unsigned int line,
377 const unsigned int face) const override
378 {
379 const static ::ndarray<unsigned int, 4, 3, 2> table = {
380 {{{{{0, 1}}, {{1, 2}}, {{2, 0}}}},
381 {{{{1, 0}}, {{0, 3}}, {{3, 1}}}},
382 {{{{0, 2}}, {{2, 3}}, {{3, 0}}}},
383 {{{{2, 1}}, {{1, 3}}, {{3, 2}}}}}};
384
385 return table[face][line];
386 }
387 };
388
389
395 {
397 vertices_of_entity(const unsigned int d,
398 const unsigned int e) const override
399 {
400 if (d == 3)
401 {
402 static const std::array<unsigned int, 5> table = {{0, 1, 2, 3, 4}};
403
404 AssertDimension(e, 0);
405
406 return {table};
407 }
408
409 if (d == 2)
410 {
411 if (e == 0)
412 {
413 static const std::array<unsigned int, 4> table = {{0, 1, 2, 3}};
414 return {table};
415 }
416
417 static const ::ndarray<unsigned int, 4, 3> table = {
418 {{{0, 2, 4}}, {{3, 1, 4}}, {{1, 0, 4}}, {{2, 3, 4}}}};
419
420 return {table[e - 1]};
421 }
422
423 if (d == 1)
424 {
425 static const ::ndarray<unsigned int, 8, 2> table = {
426 {{{0, 2}},
427 {{1, 3}},
428 {{0, 1}},
429 {{2, 3}},
430 {{0, 4}},
431 {{1, 4}},
432 {{2, 4}},
433 {{3, 4}}}};
434
435 return {table[e]};
436 }
437
439
440 return {};
441 }
442
443 virtual ReferenceCell
444 type_of_entity(const unsigned int d, const unsigned int e) const override
445 {
446 (void)e;
447
448 if (d == 3)
450
451 if (d == 2 && e == 0)
453 else if (d == 2)
455
456 if (d == 1)
458
460
462 }
463
464 unsigned int
465 n_entities(const unsigned int d) const override
466 {
467 static std::array<unsigned int, 4> table = {{5, 8, 5, 1}};
468 return table[d];
469 }
470
471 unsigned int
472 n_lines_of_surface(const unsigned int surface) const override
473 {
474 if (surface == 0)
475 return 4;
476
477 return 3;
478 }
479
480 unsigned int
481 nth_line_of_surface(const unsigned int line,
482 const unsigned int face) const override
483 {
484 const static ::ndarray<unsigned int, 5, 4> table = {
485 {{{0, 1, 2, 3}},
489 {{3, 7, 6, numbers::invalid_unsigned_int}}}};
490
491 return table[face][line];
492 }
493
494 const std::array<unsigned int, 2> &
495 vertices_of_nth_line_of_surface(const unsigned int line,
496 const unsigned int face) const override
497 {
498 static const unsigned int X = static_cast<unsigned int>(-1);
499
500 const static ::ndarray<unsigned int, 5, 4, 2> table = {
501 {{{{{0, 2}}, {{1, 3}}, {{0, 1}}, {{2, 3}}}},
502 {{{{0, 2}}, {{2, 4}}, {{4, 0}}, {{X, X}}}},
503 {{{{3, 1}}, {{1, 4}}, {{4, 3}}, {{X, X}}}},
504 {{{{1, 0}}, {{0, 4}}, {{4, 1}}, {{X, X}}}},
505 {{{{2, 3}}, {{3, 4}}, {{4, 2}}, {{X, X}}}}}};
506
507 return table[face][line];
508 }
509 };
510
511
512
517 {
519 vertices_of_entity(const unsigned int d,
520 const unsigned int e) const override
521 {
522 if (d == 3)
523 {
524 static const std::array<unsigned int, 6> table = {
525 {0, 1, 2, 3, 4, 5}};
526
527 AssertDimension(e, 0);
528
529 return {table};
530 }
531
532 if (d == 2)
533 {
534 if (e == 0 || e == 1)
535 {
536 static const ::ndarray<unsigned int, 2, 3> table = {
537 {{{1, 0, 2}}, {{3, 4, 5}}}};
538
539 return {table[e]};
540 }
541
542 static const ::ndarray<unsigned int, 3, 4> table = {
543 {{{0, 1, 3, 4}}, {{1, 2, 4, 5}}, {{2, 0, 5, 3}}}};
544
545 return {table[e - 2]};
546 }
547
548 if (d == 1)
549 {
550 static const ::ndarray<unsigned int, 9, 2> table = {
551 {{{0, 1}},
552 {{1, 2}},
553 {{2, 0}},
554 {{3, 4}},
555 {{4, 5}},
556 {{5, 3}},
557 {{0, 3}},
558 {{1, 4}},
559 {{2, 5}}}};
560
561 return {table[e]};
562 }
563
565
566 return {};
567 }
568
569 virtual ReferenceCell
570 type_of_entity(const unsigned int d, const unsigned int e) const override
571 {
572 (void)e;
573
574 if (d == 3)
576
577 if (d == 2 && e > 1)
579 else if (d == 2)
581
582 if (d == 1)
584
586
588 }
589
590 unsigned int
591 n_entities(const unsigned int d) const override
592 {
593 static std::array<unsigned int, 4> table = {{6, 9, 5, 1}};
594 return table[d];
595 }
596
597 unsigned int
598 n_lines_of_surface(const unsigned int surface) const override
599 {
600 if (surface > 1)
601 return 4;
602
603 return 3;
604 }
605
606 unsigned int
607 nth_line_of_surface(const unsigned int line,
608 const unsigned int face) const override
609 {
610 static const unsigned int X = static_cast<unsigned int>(-1);
611
612 const static ::ndarray<unsigned int, 5, 4> table = {
613 {{{0, 2, 1, X}},
614 {{3, 4, 5, X}},
615 {{6, 7, 0, 3}},
616 {{7, 8, 1, 4}},
617 {{8, 6, 5, 2}}}};
618
619 return table[face][line];
620 }
621
622 const std::array<unsigned int, 2> &
623 vertices_of_nth_line_of_surface(const unsigned int line,
624 const unsigned int face) const override
625 {
626 static const unsigned int X = static_cast<unsigned int>(-1);
627
628 const static ::ndarray<unsigned int, 5, 4, 2> table = {
629 {{{{{1, 0}}, {{0, 2}}, {{2, 1}}, {{X, X}}}},
630 {{{{3, 4}}, {{4, 5}}, {{5, 3}}, {{X, X}}}},
631 {{{{0, 3}}, {{1, 4}}, {{0, 1}}, {{3, 4}}}},
632 {{{{1, 4}}, {{2, 5}}, {{1, 2}}, {{4, 5}}}},
633 {{{{2, 5}}, {{0, 3}}, {{2, 0}}, {{5, 3}}}}}};
634
635 return table[face][line];
636 }
637 };
638
639
640
645 {
647 vertices_of_entity(const unsigned int d,
648 const unsigned int e) const override
649 {
650 if (d == 3)
651 {
652 static const std::array<unsigned int, 8> table = {
653 {0, 1, 2, 3, 4, 5, 6, 7}};
654
655 AssertDimension(e, 0);
656
657 return {table};
658 }
659
660 if (d == 2)
661 {
662 static const ::ndarray<unsigned int, 6, 4> table = {
663 {{{0, 2, 4, 6}},
664 {{1, 3, 5, 7}},
665 {{0, 4, 1, 5}},
666 {{2, 6, 3, 7}},
667 {{0, 1, 2, 3}},
668 {{4, 5, 6, 7}}}};
669
670 return {table[e]};
671 }
672
673 if (d == 1)
674 {
675 static const ::ndarray<unsigned int, 12, 2> table = {
676 {{{0, 2}},
677 {{1, 3}},
678 {{0, 1}},
679 {{2, 3}},
680 {{4, 6}},
681 {{5, 7}},
682 {{4, 5}},
683 {{6, 7}},
684 {{0, 4}},
685 {{1, 5}},
686 {{2, 6}},
687 {{3, 7}}}};
688
689 return {table[e]};
690 }
691
693
694 return {};
695 }
696
697 virtual ReferenceCell
698 type_of_entity(const unsigned int d, const unsigned int e) const override
699 {
700 (void)e;
701
702 if (d == 3)
704
705 if (d == 2)
707
708 if (d == 1)
710
712
714 }
715
716 unsigned int
717 n_entities(const unsigned int d) const override
718 {
719 static std::array<unsigned int, 4> table = {{8, 12, 6, 1}};
720 return table[d];
721 }
722
723 unsigned int
724 n_lines_of_surface(const unsigned int surface) const override
725 {
726 (void)surface;
727 return 4;
728 }
729
730 unsigned int
731 nth_line_of_surface(const unsigned int line,
732 const unsigned int face) const override
733 {
734 const static ::ndarray<unsigned int, 6, 4> table = {
735 {{{8, 10, 0, 4}},
736 {{9, 11, 1, 5}},
737 {{2, 6, 8, 9}},
738 {{3, 7, 10, 11}},
739 {{0, 1, 2, 3}},
740 {{4, 5, 6, 7}}}};
741
742 return table[face][line];
743 }
744
745 const std::array<unsigned int, 2> &
746 vertices_of_nth_line_of_surface(const unsigned int line,
747 const unsigned int face) const override
748 {
749 const static ::ndarray<unsigned int, 6, 4, 2> table = {
750 {{{{{0, 4}}, {{2, 6}}, {{0, 2}}, {{4, 6}}}},
751 {{{{1, 5}}, {{3, 7}}, {{1, 3}}, {{5, 7}}}},
752 {{{{0, 1}}, {{4, 5}}, {{0, 4}}, {{1, 5}}}},
753 {{{{2, 3}}, {{6, 7}}, {{2, 6}}, {{3, 7}}}},
754 {{{{0, 2}}, {{1, 3}}, {{0, 1}}, {{2, 3}}}},
755 {{{{4, 6}}, {{5, 7}}, {{4, 5}}, {{6, 7}}}}}};
756
757 return table[face][line];
758 }
759 };
760
761
762
769 template <typename T = unsigned int>
770 struct CRS
771 {
776 : ptr{0} {};
777
781 CRS(const std::vector<std::size_t> &ptr, const std::vector<T> &col)
782 : ptr(ptr)
783 , col(col)
784 {}
785
786 // row index
787 std::vector<std::size_t> ptr;
788
789 // column index
790 std::vector<T> col;
791 };
792
793
794
823 template <typename T = unsigned int>
825 {
826 Connectivity(const unsigned int dim,
827 const std::vector<ReferenceCell> &cell_types)
828 : dim(dim)
830 {}
831
833 entity_orientations(const unsigned int structdim)
834 {
835 if (structdim == 1)
836 return line_orientation;
837
838 AssertDimension(structdim, 2);
839
840 return quad_orientation;
841 }
842
843 inline const TriaObjectsOrientations &
844 entity_orientations(const unsigned int structdim) const
845 {
846 if (structdim == 1)
847 return line_orientation;
848
849 AssertDimension(structdim, 2);
850
851 return quad_orientation;
852 }
853
854 inline std::vector<ReferenceCell> &
855 entity_types(const unsigned int structdim)
856 {
857 if (structdim == dim)
858 return cell_types;
859
860 // for vertices/lines the entity types are clear (0/1)
861 AssertDimension(structdim, 2);
863
864 return quad_types;
865 }
866
867 inline const std::vector<ReferenceCell> &
868 entity_types(const unsigned int structdim) const
869 {
870 if (structdim == dim)
871 return cell_types;
872
873 // for vertices/lines the entity types are clear (0/1)
874 AssertDimension(structdim, 2);
876
877 return quad_types;
878 }
879
880 inline CRS<T> &
881 entity_to_entities(const unsigned int from, const unsigned int to)
882 {
883 if (from == dim && to == dim)
884 return neighbors;
885 else if (from == dim && to == dim - 1)
886 return cell_entities;
887 else if (dim == 3 && from == 2 && to == 0)
888 return quad_vertices;
889 else if (dim == 3 && from == 2 && to == 1)
890 return quad_lines;
891 else if (from == 1 && to == 0)
892 return line_vertices;
893
895
896 return cell_entities;
897 }
898
899 inline const CRS<T> &
900 entity_to_entities(const unsigned int from, const unsigned int to) const
901 {
902 if (from == dim && to == dim)
903 return neighbors;
904 else if (from == dim && to == dim - 1)
905 return cell_entities;
906 else if (dim == 3 && from == 2 && to == 0)
907 return quad_vertices;
908 else if (dim == 3 && from == 2 && to == 1)
909 return quad_lines;
910 else if (from == 1 && to == 0)
911 return line_vertices;
912
914
915 return cell_entities;
916 }
917
918 private:
919 const unsigned int dim;
920 std::vector<ReferenceCell> cell_types;
921
923
925
928
930
933
934 std::vector<ReferenceCell> quad_types;
935 };
936
937
938
946 template <typename T>
947 void
948 determine_neighbors(const CRS<T> &con_cf, CRS<T> &con_cc)
949 {
950 const auto &col_cf = con_cf.col;
951 const auto &ptr_cf = con_cf.ptr;
952
953 auto &col_cc = con_cc.col;
954 auto &ptr_cc = con_cc.ptr;
955
956 const unsigned int n_faces =
957 *std::max_element(col_cf.begin(), col_cf.end()) + 1;
958
959 // clear and initialize with -1 (assume that all faces are at the
960 // boundary)
961 col_cc = std::vector<T>(col_cf.size(), -1);
962 ptr_cc = ptr_cf;
963
964 std::vector<std::pair<T, unsigned int>> neighbors(n_faces, {-1, -1});
965
966 // loop over all cells
967 for (unsigned int i_0 = 0; i_0 < ptr_cf.size() - 1; ++i_0)
968 {
969 // ... and all its faces
970 for (std::size_t j_0 = ptr_cf[i_0]; j_0 < ptr_cf[i_0 + 1]; ++j_0)
971 {
972 if (neighbors[col_cf[j_0]].first == static_cast<unsigned int>(-1))
973 {
974 // face is visited the first time -> save the visiting cell
975 // and the face pointer
976 neighbors[col_cf[j_0]] = std::pair<T, unsigned int>(i_0, j_0);
977 }
978 else
979 {
980 // face is visited the second time -> now we know the cells
981 // on both sides of the face and we can determine for both
982 // cells the neighbor
983 col_cc[j_0] = neighbors[col_cf[j_0]].first;
984 col_cc[neighbors[col_cf[j_0]].second] = i_0;
985 }
986 }
987 }
988 }
989
990
991
999 template <int max_n_vertices, typename FU>
1000 void
1002 const unsigned int face_dimensionality,
1003 const std::vector<std::shared_ptr<CellTypeBase>> &cell_types,
1004 const std::vector<ReferenceCell> &cell_types_index,
1005 const CRS<unsigned int> &crs,
1006 CRS<unsigned int> &crs_d, // result
1007 CRS<unsigned int> &crs_0, // result
1008 TriaObjectsOrientations &orientations, // result
1009 const FU &second_key_function)
1010 {
1011 const bool compatibility_mode = true;
1012
1013 const std::vector<std::size_t> &cell_ptr = crs.ptr;
1014 const std::vector<unsigned int> &cell_vertices = crs.col;
1015 std::vector<std::size_t> &ptr_d = crs_d.ptr;
1016 std::vector<unsigned int> &col_d = crs_d.col;
1017
1018 // note: we do not pre-allocate memory for these arrays because it turned
1019 // out that counting unique entities is more expensive than push_back().
1020 std::vector<std::size_t> &ptr_0 = crs_0.ptr;
1021 std::vector<unsigned int> &col_0 = crs_0.col;
1022
1023 // clear
1024 ptr_0 = {};
1025 col_0 = {};
1026
1027 unsigned int n_entities = 0;
1028
1029 for (const auto &c : cell_types_index)
1030 n_entities += cell_types[c]->n_entities(face_dimensionality);
1031
1032 // step 1: store each d-dimensional entity of a cell (described by their
1033 // vertices) into a vector and create a key for them
1034 //
1035 // note: it turned out to be more efficient to have a vector of tuples
1036 // than to have two vectors (sorting becomes inefficient)
1037 std::vector<
1038 std::tuple<std::array<unsigned int, max_n_vertices>, unsigned int>>
1039 keys; // key (sorted vertices), cell-entity index
1040
1041 std::vector<std::array<unsigned int, max_n_vertices>> ad_entity_vertices;
1042 std::vector<ReferenceCell> ad_entity_types;
1043 std::vector<std::array<unsigned int, max_n_vertices>> ad_compatibility;
1044
1045 keys.reserve(n_entities);
1046 ad_entity_vertices.reserve(n_entities);
1047 ad_entity_types.reserve(n_entities);
1048 ad_compatibility.reserve(n_entities);
1049
1050 ptr_d.resize(cell_types_index.size() + 1);
1051 ptr_d[0] = 0;
1052
1053 static const unsigned int offset = 1;
1054
1055 // loop over all cells
1056 for (unsigned int c = 0, counter = 0; c < cell_types_index.size(); ++c)
1057 {
1058 const auto &cell_type = cell_types[cell_types_index[c]];
1059 ptr_d[c + 1] = ptr_d[c] + cell_type->n_entities(face_dimensionality);
1060
1061 // ... collect vertices of cell
1062 const ::ArrayView<const unsigned int> cell_vertice(
1063 cell_vertices.data() + cell_ptr[c], cell_ptr[c + 1] - cell_ptr[c]);
1064
1065 // ... loop over all its entities
1066 for (unsigned int e = 0;
1067 e < cell_type->n_entities(face_dimensionality);
1068 ++e)
1069 {
1070 // ... determine global entity vertices
1071 const auto &local_entity_vertices =
1072 cell_type->vertices_of_entity(face_dimensionality, e);
1073
1074 std::array<unsigned int, max_n_vertices> entity_vertices;
1075 std::fill(entity_vertices.begin(), entity_vertices.end(), 0);
1076
1077 for (unsigned int i = 0; i < local_entity_vertices.size(); ++i)
1078 entity_vertices[i] =
1079 cell_vertice[local_entity_vertices[i]] + offset;
1080
1081 // ... create key
1082 std::array<unsigned int, max_n_vertices> key = entity_vertices;
1083 std::sort(key.begin(), key.end());
1084 keys.emplace_back(key, counter++);
1085
1086 ad_entity_vertices.emplace_back(entity_vertices);
1087
1088 ad_entity_types.emplace_back(
1089 cell_type->type_of_entity(face_dimensionality, e));
1090
1091 if (compatibility_mode)
1092 ad_compatibility.emplace_back(
1093 second_key_function(entity_vertices, cell_type, c, e));
1094 }
1095 }
1096
1097 col_d.resize(keys.size());
1098 orientations.reinit(keys.size());
1099
1100 // step 2: sort according to key so that entities with same key can be
1101 // merged
1102 std::sort(keys.begin(), keys.end());
1103
1104
1105 if (compatibility_mode)
1106 {
1107 unsigned int n_unique_entities = 0;
1108 unsigned int n_unique_entity_vertices = 0;
1109
1110 std::array<unsigned int, max_n_vertices> ref_key, new_key;
1111 std::fill(ref_key.begin(), ref_key.end(), 0);
1112 for (unsigned int i = 0; i < keys.size(); ++i)
1113 {
1114 const auto offset_i = std::get<1>(keys[i]);
1115
1116 if (ref_key != std::get<0>(keys[i]))
1117 {
1118 ref_key = std::get<0>(keys[i]);
1119
1120 ++n_unique_entities;
1121 n_unique_entity_vertices +=
1122 cell_types[ad_entity_types[offset_i]]->n_entities(0);
1123
1124 new_key = ad_compatibility[offset_i];
1125 }
1126
1127 std::get<0>(keys[i]) = new_key;
1128 }
1129
1130 std::sort(keys.begin(), keys.end());
1131
1132 ptr_0.reserve(n_unique_entities + 1);
1133 col_0.reserve(n_unique_entity_vertices);
1134 }
1135
1136
1137 std::array<unsigned int, max_n_vertices> ref_key;
1138 std::array<unsigned int, max_n_vertices> ref_indices;
1139 std::fill(ref_key.begin(), ref_key.end(), 0);
1140
1141 unsigned int counter = ::numbers::invalid_unsigned_int;
1142 for (unsigned int i = 0; i < keys.size(); i++)
1143 {
1144 const auto offset_i = std::get<1>(keys[i]);
1145
1146 if (ref_key != std::get<0>(keys[i]))
1147 {
1148 // new key: default orientation is correct
1149 ++counter;
1150 ref_key = std::get<0>(keys[i]);
1151 ref_indices = ad_entity_vertices[offset_i];
1152
1153 ptr_0.push_back(col_0.size());
1154 for (const auto j : ad_entity_vertices[offset_i])
1155 if (j != 0)
1156 col_0.push_back(j - offset);
1157 }
1158 else
1159 {
1160 // previously seen key: set orientation relative to the first
1161 // occurrence
1162 orientations.set_combined_orientation(
1163 offset_i,
1164 ad_entity_types[offset_i]
1165 .template get_combined_orientation<unsigned int>(
1166 make_array_view(ad_entity_vertices[offset_i].begin(),
1167 ad_entity_vertices[offset_i].begin() +
1168 ad_entity_types[offset_i].n_vertices()),
1169 make_array_view(ref_indices.begin(),
1170 ref_indices.begin() +
1171 ad_entity_types[offset_i].n_vertices())));
1172 }
1173 col_d[offset_i] = counter;
1174 }
1175 ptr_0.push_back(col_0.size());
1176 }
1177
1178
1179
1184 template <typename FU>
1185 void
1187 const unsigned int face_dimensionality,
1188 const std::vector<std::shared_ptr<CellTypeBase>> &cell_types,
1189 const std::vector<ReferenceCell> &cell_types_index,
1190 const CRS<unsigned int> &crs,
1191 CRS<unsigned int> &crs_d,
1192 CRS<unsigned int> &crs_0,
1193 TriaObjectsOrientations &orientations,
1194 const FU &second_key_function)
1195 {
1196 std::size_t max_n_vertices = 0;
1197
1198 for (const auto &c : cell_types_index)
1199 {
1200 const auto &cell_type = cell_types[c];
1201 for (unsigned int e = 0;
1202 e < cell_type->n_entities(face_dimensionality);
1203 ++e)
1204 max_n_vertices = std::max(
1205 max_n_vertices,
1206 cell_type->vertices_of_entity(face_dimensionality, e).size());
1207 }
1208
1209 if (max_n_vertices == 2)
1210 build_face_entities_templated<2>(face_dimensionality,
1211 cell_types,
1212 cell_types_index,
1213 crs,
1214 crs_d,
1215 crs_0,
1216 orientations,
1217 second_key_function);
1218 else if (max_n_vertices == 3)
1219 build_face_entities_templated<3>(face_dimensionality,
1220 cell_types,
1221 cell_types_index,
1222 crs,
1223 crs_d,
1224 crs_0,
1225 orientations,
1226 second_key_function);
1227 else if (max_n_vertices == 4)
1228 build_face_entities_templated<4>(face_dimensionality,
1229 cell_types,
1230 cell_types_index,
1231 crs,
1232 crs_d,
1233 crs_0,
1234 orientations,
1235 second_key_function);
1236 else
1238 }
1239
1240
1241
1249 inline void
1251 const std::vector<std::shared_ptr<CellTypeBase>> &cell_types,
1252 const std::vector<ReferenceCell> &cell_types_index,
1253 const CRS<unsigned int> &con_cv,
1254 const CRS<unsigned int> &con_cl,
1255 const CRS<unsigned int> &con_lv,
1256 const CRS<unsigned int> &con_cq,
1257 const CRS<unsigned int> &con_qv,
1258 const TriaObjectsOrientations &ori_cq,
1259 CRS<unsigned int> &con_ql, // result
1260 TriaObjectsOrientations &ori_ql, // result
1261 std::vector<ReferenceCell> &quad_t_id // result
1262 )
1263 {
1264 // reset output
1265 con_ql.ptr = {};
1266 con_ql.col = {};
1267
1268 con_ql.ptr.resize(con_qv.ptr.size());
1269 con_ql.ptr[0] = 0;
1270
1271 quad_t_id.resize(con_qv.ptr.size() - 1);
1272
1273 // count the number of lines of each face
1274 for (unsigned int c = 0; c < con_cq.ptr.size() - 1; ++c)
1275 {
1276 const auto &cell_type = cell_types[cell_types_index[c]];
1277
1278 // loop over faces
1279 for (unsigned int f_ = con_cq.ptr[c], f_index = 0;
1280 f_ < con_cq.ptr[c + 1];
1281 ++f_, ++f_index)
1282 {
1283 const unsigned int f = con_cq.col[f_];
1284
1285 con_ql.ptr[f + 1] = cell_type->n_lines_of_surface(f_index);
1286 }
1287 }
1288
1289 // use the counts to determine the offsets -> prefix sum
1290 for (unsigned int i = 0; i < con_ql.ptr.size() - 1; ++i)
1291 con_ql.ptr[i + 1] += con_ql.ptr[i];
1292
1293 // allocate memory
1294 con_ql.col.resize(con_ql.ptr.back());
1295 ori_ql.reinit(con_ql.ptr.back());
1296
1297 // loop over cells
1298 for (unsigned int c = 0; c < con_cq.ptr.size() - 1; ++c)
1299 {
1300 const auto &cell_type = cell_types[cell_types_index[c]];
1301
1302 // loop over faces
1303 for (unsigned int f_ = con_cq.ptr[c], f_index = 0;
1304 f_ < con_cq.ptr[c + 1];
1305 ++f_, ++f_index)
1306 {
1307 const unsigned int f = con_cq.col[f_];
1308
1309 // only faces with default orientation have to do something
1310 if (ori_cq.get_combined_orientation(f_) !=
1312 continue;
1313
1314 // determine entity type of face
1315 quad_t_id[f] = cell_type->type_of_entity(2, f_index);
1316
1317 // loop over lines
1318 for (unsigned int l = 0;
1319 l < cell_type->n_lines_of_surface(f_index);
1320 ++l)
1321 {
1322 // determine global index of line
1323 const unsigned int local_line_index =
1324 cell_type->nth_line_of_surface(l, f_index);
1325 const unsigned int global_line_index =
1326 con_cl.col[con_cl.ptr[c] + local_line_index];
1327 con_ql.col[con_ql.ptr[f] + l] = global_line_index;
1328
1329 // determine orientation of line
1330 const auto line_vertices_1_ref =
1331 cell_type->vertices_of_nth_line_of_surface(l, f_index);
1332
1333 bool same = true;
1334 for (unsigned int v = 0; v < line_vertices_1_ref.size(); ++v)
1335 if (con_cv.col[con_cv.ptr[c] + line_vertices_1_ref[v]] !=
1336 con_lv.col[con_lv.ptr[global_line_index] + v])
1337 {
1338 same = false;
1339 break;
1340 }
1341
1342 // ... comparison gives orientation
1344 con_ql.ptr[f] + l,
1347 }
1348 }
1349 }
1350 }
1351
1352
1353
1362 template <typename T>
1363 Connectivity<T>
1364 build_connectivity(const unsigned int dim,
1365 const std::vector<std::shared_ptr<CellTypeBase>> &cell_t,
1366 const std::vector<ReferenceCell> &cell_t_id,
1367 const CRS<T> &con_cv)
1368 {
1369 Connectivity<T> connectivity(dim, cell_t_id);
1370
1371 CRS<T> temp1; // needed for 3d
1372
1373 if (dim == 1)
1374 connectivity.entity_to_entities(1, 0) = con_cv;
1375
1376 if (dim == 2 || dim == 3) // build lines
1377 {
1379
1381 1,
1382 cell_t,
1383 connectivity.entity_types(dim),
1384 con_cv,
1385 dim == 2 ? connectivity.entity_to_entities(2, 1) : temp1,
1386 connectivity.entity_to_entities(1, 0),
1387 dim == 2 ? connectivity.entity_orientations(1) : dummy,
1388 [](auto key, const auto &, const auto &, const auto &) {
1389 // to ensure same enumeration as in deal.II
1390 return key;
1391 });
1392 }
1393
1394 if (dim == 3) // build quads
1395 {
1397 2,
1398 cell_t,
1399 connectivity.entity_types(3),
1400 con_cv,
1401 connectivity.entity_to_entities(3, 2),
1402 connectivity.entity_to_entities(2, 0),
1403 connectivity.entity_orientations(2),
1404 [&](auto key, const auto &cell_type, const auto &c, const auto &f) {
1405 // to ensure same enumeration as in deal.II
1406 AssertIndexRange(cell_type->n_lines_of_surface(f),
1407 key.size() + 1);
1408
1409 unsigned int l = 0;
1410
1411 for (; l < cell_type->n_lines_of_surface(f); ++l)
1412 key[l] =
1413 temp1
1414 .col[temp1.ptr[c] + cell_type->nth_line_of_surface(l, f)] +
1415 1 /*offset!*/;
1416
1417 for (; l < key.size(); ++l)
1418 key[l] = 0;
1419
1420 return key;
1421 });
1422
1423 // create connectivity: quad -> line
1424 build_intersection(cell_t,
1425 connectivity.entity_types(3),
1426 con_cv,
1427 temp1,
1428 connectivity.entity_to_entities(1, 0),
1429 connectivity.entity_to_entities(3, 2),
1430 connectivity.entity_to_entities(2, 0),
1431 connectivity.entity_orientations(2),
1432 connectivity.entity_to_entities(2, 1),
1433 connectivity.entity_orientations(1),
1434 connectivity.entity_types(2));
1435 }
1436
1437 // determine neighbors
1438 determine_neighbors(connectivity.entity_to_entities(dim, dim - 1),
1439 connectivity.entity_to_entities(dim, dim));
1440
1441 return connectivity;
1442 }
1443
1444
1445
1449 template <typename T, int dim>
1450 Connectivity<T>
1451 build_connectivity(const std::vector<CellData<dim>> &cells)
1452 {
1453 AssertThrow(cells.size() > 0, ExcMessage("No cells have been provided!"));
1454
1455 // vector of possible cell entity types
1456 std::vector<std::shared_ptr<CellTypeBase>> cell_types_impl(8);
1457
1458 cell_types_impl[ReferenceCells::Line] = std::make_shared<CellTypeLine>();
1459 cell_types_impl[ReferenceCells::Triangle] =
1460 std::make_shared<CellTypeTriangle>();
1461 cell_types_impl[ReferenceCells::Quadrilateral] =
1462 std::make_shared<CellTypeQuadrilateral>();
1463 cell_types_impl[ReferenceCells::Tetrahedron] =
1464 std::make_shared<CellTypeTetrahedron>();
1465 cell_types_impl[ReferenceCells::Pyramid] =
1466 std::make_shared<CellTypePyramid>();
1467 cell_types_impl[ReferenceCells::Wedge] =
1468 std::make_shared<CellTypeWedge>();
1469 cell_types_impl[ReferenceCells::Hexahedron] =
1470 std::make_shared<CellTypeHexahedron>();
1471
1472 // determine cell types and process vertices
1473 std::vector<T> cell_vertices;
1474 cell_vertices.reserve(
1475 std::accumulate(cells.begin(),
1476 cells.end(),
1477 0,
1478 [](const auto &result, const auto &cell) {
1479 return result + cell.vertices.size();
1480 }));
1481
1482 std::vector<std::size_t> cell_vertices_ptr;
1483 cell_vertices_ptr.reserve(cells.size() + 1);
1484 cell_vertices_ptr.push_back(0);
1485
1486 std::vector<ReferenceCell> cell_types_indices;
1487 cell_types_indices.reserve(cells.size());
1488
1489 // loop over cells and create CRS
1490 for (const auto &cell : cells)
1491 {
1492#ifdef DEBUG
1493 auto vertices_unique = cell.vertices;
1494 std::sort(vertices_unique.begin(), vertices_unique.end());
1495 vertices_unique.erase(std::unique(vertices_unique.begin(),
1496 vertices_unique.end()),
1497 vertices_unique.end());
1498
1499 Assert(vertices_unique.size() == cell.vertices.size(),
1500 ExcMessage(
1501 "The definition of a cell refers to the same vertex several "
1502 "times. This is not possible. A common reason is that "
1503 "CellData::vertices has a size that does not match the "
1504 "size expected from the reference cell. Please resize "
1505 "CellData::vertices or use the appropriate constructor of "
1506 "CellData."));
1507#endif
1508
1509 const ReferenceCell reference_cell =
1510 ReferenceCell::n_vertices_to_type(dim, cell.vertices.size());
1511
1512 Assert(reference_cell != ReferenceCells::Invalid,
1514 AssertIndexRange(static_cast<std::uint8_t>(reference_cell),
1515 cell_types_impl.size());
1516 Assert(cell_types_impl[reference_cell].get() != nullptr,
1518
1519 cell_types_indices.push_back(reference_cell);
1520
1521 // create CRS of vertices (to remove template argument dim)
1522 for (const auto &vertex : cell.vertices)
1523 cell_vertices.push_back(vertex);
1524
1525 cell_vertices_ptr.push_back(cell_vertices.size());
1526 }
1527
1528 // do the actual work
1529 return build_connectivity<T>(dim,
1530 cell_types_impl,
1531 cell_types_indices,
1532 {cell_vertices_ptr, cell_vertices});
1533 }
1534 } // namespace TriangulationImplementation
1535} // namespace internal
1536
1537
1539
1540#endif
ArrayView< std::remove_reference_t< typename std::iterator_traits< Iterator >::reference >, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition array_view.h:949
static ReferenceCell n_vertices_to_type(const int dim, const unsigned int n_vertices)
types::geometric_orientation get_combined_orientation(const unsigned int object) const
void set_combined_orientation(const unsigned int object, const types::geometric_orientation value)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:518
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:519
Point< 2 > first
Definition grid_out.cc:4629
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
#define DEAL_II_NOT_IMPLEMENTED()
constexpr const ReferenceCell Tetrahedron
constexpr const ReferenceCell Quadrilateral
constexpr const ReferenceCell Wedge
constexpr const ReferenceCell Pyramid
constexpr const ReferenceCell Invalid
constexpr const ReferenceCell Triangle
constexpr const ReferenceCell Hexahedron
constexpr const ReferenceCell Vertex
constexpr const ReferenceCell Line
void determine_neighbors(const CRS< T > &con_cf, CRS< T > &con_cc)
void build_intersection(const std::vector< std::shared_ptr< CellTypeBase > > &cell_types, const std::vector< ReferenceCell > &cell_types_index, const CRS< unsigned int > &con_cv, const CRS< unsigned int > &con_cl, const CRS< unsigned int > &con_lv, const CRS< unsigned int > &con_cq, const CRS< unsigned int > &con_qv, const TriaObjectsOrientations &ori_cq, CRS< unsigned int > &con_ql, TriaObjectsOrientations &ori_ql, std::vector< ReferenceCell > &quad_t_id)
void build_face_entities_templated(const unsigned int face_dimensionality, const std::vector< std::shared_ptr< CellTypeBase > > &cell_types, const std::vector< ReferenceCell > &cell_types_index, const CRS< unsigned int > &crs, CRS< unsigned int > &crs_d, CRS< unsigned int > &crs_0, TriaObjectsOrientations &orientations, const FU &second_key_function)
Connectivity< T > build_connectivity(const unsigned int dim, const std::vector< std::shared_ptr< CellTypeBase > > &cell_t, const std::vector< ReferenceCell > &cell_t_id, const CRS< T > &con_cv)
void build_face_entities(const unsigned int face_dimensionality, const std::vector< std::shared_ptr< CellTypeBase > > &cell_types, const std::vector< ReferenceCell > &cell_types_index, const CRS< unsigned int > &crs, CRS< unsigned int > &crs_d, CRS< unsigned int > &crs_0, TriaObjectsOrientations &orientations, const FU &second_key_function)
constexpr unsigned int invalid_unsigned_int
Definition types.h:232
constexpr types::geometric_orientation reverse_line_orientation
Definition types.h:359
constexpr types::geometric_orientation default_geometric_orientation
Definition types.h:346
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
CRS(const std::vector< std::size_t > &ptr, const std::vector< T > &col)
virtual unsigned int n_entities(const unsigned int d) const
virtual ::ArrayView< const unsigned int > vertices_of_entity(const unsigned int d, const unsigned int e) const
virtual ReferenceCell type_of_entity(const unsigned int d, const unsigned int e) const
virtual const std::array< unsigned int, 2 > & vertices_of_nth_line_of_surface(const unsigned int line, const unsigned int face) const
virtual unsigned int n_lines_of_surface(const unsigned int face) const
virtual unsigned int nth_line_of_surface(const unsigned int line, const unsigned int face) const
virtual ReferenceCell type_of_entity(const unsigned int d, const unsigned int e) const override
::ArrayView< const unsigned int > vertices_of_entity(const unsigned int d, const unsigned int e) const override
unsigned int nth_line_of_surface(const unsigned int line, const unsigned int face) const override
unsigned int n_entities(const unsigned int d) const override
const std::array< unsigned int, 2 > & vertices_of_nth_line_of_surface(const unsigned int line, const unsigned int face) const override
unsigned int n_lines_of_surface(const unsigned int surface) const override
::ArrayView< const unsigned int > vertices_of_entity(const unsigned int d, const unsigned int e) const override
unsigned int n_entities(const unsigned int d) const override
ReferenceCell type_of_entity(const unsigned int d, const unsigned int e) const override
unsigned int n_lines_of_surface(const unsigned int surface) const override
const std::array< unsigned int, 2 > & vertices_of_nth_line_of_surface(const unsigned int line, const unsigned int face) const override
::ArrayView< const unsigned int > vertices_of_entity(const unsigned int d, const unsigned int e) const override
unsigned int n_entities(const unsigned int d) const override
unsigned int nth_line_of_surface(const unsigned int line, const unsigned int face) const override
virtual ReferenceCell type_of_entity(const unsigned int d, const unsigned int e) const override
virtual ReferenceCell type_of_entity(const unsigned int d, const unsigned int e) const override
unsigned int n_entities(const unsigned int d) const override
::ArrayView< const unsigned int > vertices_of_entity(const unsigned int d, const unsigned int e) const override
unsigned int nth_line_of_surface(const unsigned int line, const unsigned int face) const override
unsigned int n_lines_of_surface(const unsigned int line) const override
::ArrayView< const unsigned int > vertices_of_entity(const unsigned int d, const unsigned int e) const override
unsigned int n_entities(const unsigned int d) const override
const std::array< unsigned int, 2 > & vertices_of_nth_line_of_surface(const unsigned int line, const unsigned int face) const override
virtual ReferenceCell type_of_entity(const unsigned int d, const unsigned int e) const override
unsigned int n_entities(const unsigned int d) const override
::ArrayView< const unsigned int > vertices_of_entity(const unsigned int d, const unsigned int e) const override
virtual ReferenceCell type_of_entity(const unsigned int d, const unsigned int e) const override
unsigned int n_entities(const unsigned int d) const override
unsigned int n_lines_of_surface(const unsigned int surface) const override
::ArrayView< const unsigned int > vertices_of_entity(const unsigned int d, const unsigned int e) const override
virtual ReferenceCell type_of_entity(const unsigned int d, const unsigned int e) const override
unsigned int nth_line_of_surface(const unsigned int line, const unsigned int face) const override
const std::array< unsigned int, 2 > & vertices_of_nth_line_of_surface(const unsigned int line, const unsigned int face) const override
const CRS< T > & entity_to_entities(const unsigned int from, const unsigned int to) const
CRS< T > & entity_to_entities(const unsigned int from, const unsigned int to)
std::vector< ReferenceCell > & entity_types(const unsigned int structdim)
const std::vector< ReferenceCell > & entity_types(const unsigned int structdim) const
Connectivity(const unsigned int dim, const std::vector< ReferenceCell > &cell_types)
TriaObjectsOrientations & entity_orientations(const unsigned int structdim)
const TriaObjectsOrientations & entity_orientations(const unsigned int structdim) const