Reference documentation for deal.II version GIT relicensing-464-g14f7274e4d 2024-04-22 15:40:02+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
29
30
31namespace internal
32{
33 namespace TriangulationImplementation
34 {
40 {
44 virtual ~CellTypeBase() = default;
45
49 virtual unsigned int
50 n_entities(const unsigned int d) const
51 {
53 (void)d;
54
55 return 0;
56 }
57
61 virtual ::ArrayView<const unsigned int>
62 vertices_of_entity(const unsigned int d, const unsigned int e) const
63 {
65 (void)d;
66 (void)e;
67
68 return {};
69 }
70
74 virtual ReferenceCell
75 type_of_entity(const unsigned int d, const unsigned int e) const
76 {
78 (void)d;
79 (void)e;
80
82 }
83
87 virtual unsigned int
88 n_lines_of_surface(const unsigned int face) const
89 {
91 (void)face;
92
93 return 0;
94 }
95
99 virtual unsigned int
100 nth_line_of_surface(const unsigned int line,
101 const unsigned int face) const
102 {
104 (void)line;
105 (void)face;
106
107 return 0;
108 }
109
113 virtual const std::array<unsigned int, 2> &
114 vertices_of_nth_line_of_surface(const unsigned int line,
115 const unsigned int face) const
116 {
118 (void)line;
119 (void)face;
120
121 const static std::array<unsigned int, 2> table = {};
122
123 return table;
124 }
125 };
126
127
128
133 {
135 vertices_of_entity(const unsigned int d,
136 const unsigned int e) const override
137 {
138 (void)e;
139
140 if (d == 1)
141 {
142 static const std::array<unsigned int, 2> table = {{0, 1}};
143
144 AssertDimension(e, 0);
145
146 return {table};
147 }
148
150
151 return {};
152 }
153
155 type_of_entity(const unsigned int d, const unsigned int e) const override
156 {
157 (void)e;
158
159 if (d == 1)
161
163
165 }
166
167 unsigned int
168 n_entities(const unsigned int d) const override
169 {
170 static std::array<unsigned int, 2> table = {{2, 1}};
171 return table[d];
172 }
173 };
174
175
176
181 {
183 vertices_of_entity(const unsigned int d,
184 const unsigned int e) const override
185 {
186 if (d == 2)
187 {
188 static const std::array<unsigned int, 3> table = {{0, 1, 2}};
189
190 AssertDimension(e, 0);
191
192 return {table};
193 }
194
195 if (d == 1)
196 {
197 static const ::ndarray<unsigned int, 3, 2> table = {
198 {{{0, 1}}, {{1, 2}}, {{2, 0}}}};
199
200 return {table[e]};
201 }
202
204
205 return {};
206 }
207
208 virtual ReferenceCell
209 type_of_entity(const unsigned int d, const unsigned int e) const override
210 {
211 (void)e;
212
213 if (d == 2)
215
216 if (d == 1)
218
220
222 }
223
224 unsigned int
225 n_entities(const unsigned int d) const override
226 {
227 static std::array<unsigned int, 3> table = {{3, 3, 1}};
228 return table[d];
229 }
230 };
231
232
233
238 {
240 vertices_of_entity(const unsigned int d,
241 const unsigned int e) const override
242 {
243 if (d == 2)
244 {
245 static const std::array<unsigned int, 4> table = {{0, 1, 2, 3}};
246
247 AssertDimension(e, 0);
248
249 return {table};
250 }
251
252 if (d == 1)
253 {
254 static const ::ndarray<unsigned int, 4, 2> table = {
255 {{{0, 2}}, {{1, 3}}, {{0, 1}}, {{2, 3}}}};
256
257 return {table[e]};
258 }
259
261
262 return {};
263 }
264
265 virtual ReferenceCell
266 type_of_entity(const unsigned int d, const unsigned int e) const override
267 {
268 (void)e;
269
270 if (d == 2)
272
273 if (d == 1)
275
277
279 }
280
281 unsigned int
282 n_entities(const unsigned int d) const override
283 {
284 static std::array<unsigned int, 3> table = {{4, 4, 1}};
285 return table[d];
286 }
287 };
288
289
290
295 {
297 vertices_of_entity(const unsigned int d,
298 const unsigned int e) const override
299 {
300 if (d == 3)
301 {
302 static const std::array<unsigned int, 4> table = {{0, 1, 2, 3}};
303
304 AssertDimension(e, 0);
305
306 return {table};
307 }
308
309 if (d == 2)
310 {
311 static const ::ndarray<unsigned int, 4, 3> table = {
312 {{{0, 1, 2}}, {{1, 0, 3}}, {{0, 2, 3}}, {{2, 1, 3}}}};
313
314 return {table[e]};
315 }
316
317 if (d == 1)
318 {
319 static const ::ndarray<unsigned int, 6, 2> table = {
320 {{{0, 1}}, {{1, 2}}, {{2, 0}}, {{0, 3}}, {{1, 3}}, {{2, 3}}}};
321
322 return {table[e]};
323 }
324
326
327 return {};
328 }
329
330 virtual ReferenceCell
331 type_of_entity(const unsigned int d, const unsigned int e) const override
332 {
333 (void)e;
334
335 if (d == 3)
337
338 if (d == 2)
340
341 if (d == 1)
343
345
347 }
348
349 unsigned int
350 n_entities(const unsigned int d) const override
351 {
352 static std::array<unsigned int, 4> table = {{4, 6, 4, 1}};
353 return table[d];
354 }
355
356 unsigned int
357 n_lines_of_surface(const unsigned int line) const override
358 {
359 (void)line;
360 return 3;
361 }
362
363 unsigned int
364 nth_line_of_surface(const unsigned int line,
365 const unsigned int face) const override
366 {
367 const static ::ndarray<unsigned int, 4, 3> table = {
368 {{{0, 1, 2}}, {{0, 3, 4}}, {{2, 5, 3}}, {{1, 4, 5}}}};
369
370 return table[face][line];
371 }
372
373 const std::array<unsigned int, 2> &
374 vertices_of_nth_line_of_surface(const unsigned int line,
375 const unsigned int face) const override
376 {
377 const static ::ndarray<unsigned int, 4, 3, 2> table = {
378 {{{{{0, 1}}, {{1, 2}}, {{2, 0}}}},
379 {{{{1, 0}}, {{0, 3}}, {{3, 1}}}},
380 {{{{0, 2}}, {{2, 3}}, {{3, 0}}}},
381 {{{{2, 1}}, {{1, 3}}, {{3, 2}}}}}};
382
383 return table[face][line];
384 }
385 };
386
387
393 {
395 vertices_of_entity(const unsigned int d,
396 const unsigned int e) const override
397 {
398 if (d == 3)
399 {
400 static const std::array<unsigned int, 5> table = {{0, 1, 2, 3, 4}};
401
402 AssertDimension(e, 0);
403
404 return {table};
405 }
406
407 if (d == 2)
408 {
409 if (e == 0)
410 {
411 static const std::array<unsigned int, 4> table = {{0, 1, 2, 3}};
412 return {table};
413 }
414
415 static const ::ndarray<unsigned int, 4, 3> table = {
416 {{{0, 2, 4}}, {{3, 1, 4}}, {{1, 0, 4}}, {{2, 3, 4}}}};
417
418 return {table[e - 1]};
419 }
420
421 if (d == 1)
422 {
423 static const ::ndarray<unsigned int, 8, 2> table = {
424 {{{0, 2}},
425 {{1, 3}},
426 {{0, 1}},
427 {{2, 3}},
428 {{0, 4}},
429 {{1, 4}},
430 {{2, 4}},
431 {{3, 4}}}};
432
433 return {table[e]};
434 }
435
437
438 return {};
439 }
440
441 virtual ReferenceCell
442 type_of_entity(const unsigned int d, const unsigned int e) const override
443 {
444 (void)e;
445
446 if (d == 3)
448
449 if (d == 2 && e == 0)
451 else if (d == 2)
453
454 if (d == 1)
456
458
460 }
461
462 unsigned int
463 n_entities(const unsigned int d) const override
464 {
465 static std::array<unsigned int, 4> table = {{5, 8, 5, 1}};
466 return table[d];
467 }
468
469 unsigned int
470 n_lines_of_surface(const unsigned int surface) const override
471 {
472 if (surface == 0)
473 return 4;
474
475 return 3;
476 }
477
478 unsigned int
479 nth_line_of_surface(const unsigned int line,
480 const unsigned int face) const override
481 {
482 const static ::ndarray<unsigned int, 5, 4> table = {
483 {{{0, 1, 2, 3}},
487 {{3, 7, 6, numbers::invalid_unsigned_int}}}};
488
489 return table[face][line];
490 }
491
492 const std::array<unsigned int, 2> &
493 vertices_of_nth_line_of_surface(const unsigned int line,
494 const unsigned int face) const override
495 {
496 static const unsigned int X = static_cast<unsigned int>(-1);
497
498 const static ::ndarray<unsigned int, 5, 4, 2> table = {
499 {{{{{0, 2}}, {{1, 3}}, {{0, 1}}, {{2, 3}}}},
500 {{{{0, 2}}, {{2, 4}}, {{4, 0}}, {{X, X}}}},
501 {{{{3, 1}}, {{1, 4}}, {{4, 3}}, {{X, X}}}},
502 {{{{1, 0}}, {{0, 4}}, {{4, 1}}, {{X, X}}}},
503 {{{{2, 3}}, {{3, 4}}, {{4, 2}}, {{X, X}}}}}};
504
505 return table[face][line];
506 }
507 };
508
509
510
515 {
517 vertices_of_entity(const unsigned int d,
518 const unsigned int e) const override
519 {
520 if (d == 3)
521 {
522 static const std::array<unsigned int, 6> table = {
523 {0, 1, 2, 3, 4, 5}};
524
525 AssertDimension(e, 0);
526
527 return {table};
528 }
529
530 if (d == 2)
531 {
532 if (e == 0 || e == 1)
533 {
534 static const ::ndarray<unsigned int, 2, 3> table = {
535 {{{1, 0, 2}}, {{3, 4, 5}}}};
536
537 return {table[e]};
538 }
539
540 static const ::ndarray<unsigned int, 3, 4> table = {
541 {{{0, 1, 3, 4}}, {{1, 2, 4, 5}}, {{2, 0, 5, 3}}}};
542
543 return {table[e - 2]};
544 }
545
546 if (d == 1)
547 {
548 static const ::ndarray<unsigned int, 9, 2> table = {
549 {{{0, 1}},
550 {{1, 2}},
551 {{2, 0}},
552 {{3, 4}},
553 {{4, 5}},
554 {{5, 3}},
555 {{0, 3}},
556 {{1, 4}},
557 {{2, 5}}}};
558
559 return {table[e]};
560 }
561
563
564 return {};
565 }
566
567 virtual ReferenceCell
568 type_of_entity(const unsigned int d, const unsigned int e) const override
569 {
570 (void)e;
571
572 if (d == 3)
574
575 if (d == 2 && e > 1)
577 else if (d == 2)
579
580 if (d == 1)
582
584
586 }
587
588 unsigned int
589 n_entities(const unsigned int d) const override
590 {
591 static std::array<unsigned int, 4> table = {{6, 9, 5, 1}};
592 return table[d];
593 }
594
595 unsigned int
596 n_lines_of_surface(const unsigned int surface) const override
597 {
598 if (surface > 1)
599 return 4;
600
601 return 3;
602 }
603
604 unsigned int
605 nth_line_of_surface(const unsigned int line,
606 const unsigned int face) const override
607 {
608 static const unsigned int X = static_cast<unsigned int>(-1);
609
610 const static ::ndarray<unsigned int, 5, 4> table = {
611 {{{0, 2, 1, X}},
612 {{3, 4, 5, X}},
613 {{6, 7, 0, 3}},
614 {{7, 8, 1, 4}},
615 {{8, 6, 5, 2}}}};
616
617 return table[face][line];
618 }
619
620 const std::array<unsigned int, 2> &
621 vertices_of_nth_line_of_surface(const unsigned int line,
622 const unsigned int face) const override
623 {
624 static const unsigned int X = static_cast<unsigned int>(-1);
625
626 const static ::ndarray<unsigned int, 5, 4, 2> table = {
627 {{{{{1, 0}}, {{0, 2}}, {{2, 1}}, {{X, X}}}},
628 {{{{3, 4}}, {{4, 5}}, {{5, 3}}, {{X, X}}}},
629 {{{{0, 3}}, {{1, 4}}, {{0, 1}}, {{3, 4}}}},
630 {{{{1, 4}}, {{2, 5}}, {{1, 2}}, {{4, 5}}}},
631 {{{{2, 5}}, {{0, 3}}, {{2, 0}}, {{5, 3}}}}}};
632
633 return table[face][line];
634 }
635 };
636
637
638
643 {
645 vertices_of_entity(const unsigned int d,
646 const unsigned int e) const override
647 {
648 if (d == 3)
649 {
650 static const std::array<unsigned int, 8> table = {
651 {0, 1, 2, 3, 4, 5, 6, 7}};
652
653 AssertDimension(e, 0);
654
655 return {table};
656 }
657
658 if (d == 2)
659 {
660 static const ::ndarray<unsigned int, 6, 4> table = {
661 {{{0, 2, 4, 6}},
662 {{1, 3, 5, 7}},
663 {{0, 4, 1, 5}},
664 {{2, 6, 3, 7}},
665 {{0, 1, 2, 3}},
666 {{4, 5, 6, 7}}}};
667
668 return {table[e]};
669 }
670
671 if (d == 1)
672 {
673 static const ::ndarray<unsigned int, 12, 2> table = {
674 {{{0, 2}},
675 {{1, 3}},
676 {{0, 1}},
677 {{2, 3}},
678 {{4, 6}},
679 {{5, 7}},
680 {{4, 5}},
681 {{6, 7}},
682 {{0, 4}},
683 {{1, 5}},
684 {{2, 6}},
685 {{3, 7}}}};
686
687 return {table[e]};
688 }
689
691
692 return {};
693 }
694
695 virtual ReferenceCell
696 type_of_entity(const unsigned int d, const unsigned int e) const override
697 {
698 (void)e;
699
700 if (d == 3)
702
703 if (d == 2)
705
706 if (d == 1)
708
710
712 }
713
714 unsigned int
715 n_entities(const unsigned int d) const override
716 {
717 static std::array<unsigned int, 4> table = {{8, 12, 6, 1}};
718 return table[d];
719 }
720
721 unsigned int
722 n_lines_of_surface(const unsigned int surface) const override
723 {
724 (void)surface;
725 return 4;
726 }
727
728 unsigned int
729 nth_line_of_surface(const unsigned int line,
730 const unsigned int face) const override
731 {
732 const static ::ndarray<unsigned int, 6, 4> table = {
733 {{{8, 10, 0, 4}},
734 {{9, 11, 1, 5}},
735 {{2, 6, 8, 9}},
736 {{3, 7, 10, 11}},
737 {{0, 1, 2, 3}},
738 {{4, 5, 6, 7}}}};
739
740 return table[face][line];
741 }
742
743 const std::array<unsigned int, 2> &
744 vertices_of_nth_line_of_surface(const unsigned int line,
745 const unsigned int face) const override
746 {
747 const static ::ndarray<unsigned int, 6, 4, 2> table = {
748 {{{{{0, 4}}, {{2, 6}}, {{0, 2}}, {{4, 6}}}},
749 {{{{1, 5}}, {{3, 7}}, {{1, 3}}, {{5, 7}}}},
750 {{{{0, 1}}, {{4, 5}}, {{0, 4}}, {{1, 5}}}},
751 {{{{2, 3}}, {{6, 7}}, {{2, 6}}, {{3, 7}}}},
752 {{{{0, 2}}, {{1, 3}}, {{0, 1}}, {{2, 3}}}},
753 {{{{4, 6}}, {{5, 7}}, {{4, 5}}, {{6, 7}}}}}};
754
755 return table[face][line];
756 }
757 };
758
759
760
767 template <typename T = unsigned int>
768 struct CRS
769 {
774 : ptr{0} {};
775
779 CRS(const std::vector<std::size_t> &ptr, const std::vector<T> &col)
780 : ptr(ptr)
781 , col(col)
782 {}
783
784 // row index
785 std::vector<std::size_t> ptr;
786
787 // column index
788 std::vector<T> col;
789 };
790
791
792
821 template <typename T = unsigned int>
823 {
824 Connectivity(const unsigned int dim,
825 const std::vector<ReferenceCell> &cell_types)
826 : dim(dim)
828 {}
829
831 entity_orientations(const unsigned int structdim)
832 {
833 if (structdim == 1)
834 return line_orientation;
835
836 AssertDimension(structdim, 2);
837
838 return quad_orientation;
839 }
840
841 inline const TriaObjectsOrientations &
842 entity_orientations(const unsigned int structdim) const
843 {
844 if (structdim == 1)
845 return line_orientation;
846
847 AssertDimension(structdim, 2);
848
849 return quad_orientation;
850 }
851
852 inline std::vector<ReferenceCell> &
853 entity_types(const unsigned int structdim)
854 {
855 if (structdim == dim)
856 return cell_types;
857
858 // for vertices/lines the entity types are clear (0/1)
859 AssertDimension(structdim, 2);
861
862 return quad_types;
863 }
864
865 inline const std::vector<ReferenceCell> &
866 entity_types(const unsigned int structdim) const
867 {
868 if (structdim == dim)
869 return cell_types;
870
871 // for vertices/lines the entity types are clear (0/1)
872 AssertDimension(structdim, 2);
874
875 return quad_types;
876 }
877
878 inline CRS<T> &
879 entity_to_entities(const unsigned int from, const unsigned int to)
880 {
881 if (from == dim && to == dim)
882 return neighbors;
883 else if (from == dim && to == dim - 1)
884 return cell_entities;
885 else if (dim == 3 && from == 2 && to == 0)
886 return quad_vertices;
887 else if (dim == 3 && from == 2 && to == 1)
888 return quad_lines;
889 else if (from == 1 && to == 0)
890 return line_vertices;
891
893
894 return cell_entities;
895 }
896
897 inline const CRS<T> &
898 entity_to_entities(const unsigned int from, const unsigned int to) const
899 {
900 if (from == dim && to == dim)
901 return neighbors;
902 else if (from == dim && to == dim - 1)
903 return cell_entities;
904 else if (dim == 3 && from == 2 && to == 0)
905 return quad_vertices;
906 else if (dim == 3 && from == 2 && to == 1)
907 return quad_lines;
908 else if (from == 1 && to == 0)
909 return line_vertices;
910
912
913 return cell_entities;
914 }
915
916 private:
917 const unsigned int dim;
918 std::vector<ReferenceCell> cell_types;
919
921
923
926
928
931
932 std::vector<ReferenceCell> quad_types;
933 };
934
935
936
944 template <typename T>
945 void
946 determine_neighbors(const CRS<T> &con_cf, CRS<T> &con_cc)
947 {
948 const auto &col_cf = con_cf.col;
949 const auto &ptr_cf = con_cf.ptr;
950
951 auto &col_cc = con_cc.col;
952 auto &ptr_cc = con_cc.ptr;
953
954 const unsigned int n_faces =
955 *std::max_element(col_cf.begin(), col_cf.end()) + 1;
956
957 // clear and initialize with -1 (assume that all faces are at the
958 // boundary)
959 col_cc = std::vector<T>(col_cf.size(), -1);
960 ptr_cc = ptr_cf;
961
962 std::vector<std::pair<T, unsigned int>> neighbors(n_faces, {-1, -1});
963
964 // loop over all cells
965 for (unsigned int i_0 = 0; i_0 < ptr_cf.size() - 1; ++i_0)
966 {
967 // ... and all its faces
968 for (std::size_t j_0 = ptr_cf[i_0]; j_0 < ptr_cf[i_0 + 1]; ++j_0)
969 {
970 if (neighbors[col_cf[j_0]].first == static_cast<unsigned int>(-1))
971 {
972 // face is visited the first time -> save the visiting cell
973 // and the face pointer
974 neighbors[col_cf[j_0]] = std::pair<T, unsigned int>(i_0, j_0);
975 }
976 else
977 {
978 // face is visited the second time -> now we know the cells
979 // on both sides of the face and we can determine for both
980 // cells the neighbor
981 col_cc[j_0] = neighbors[col_cf[j_0]].first;
982 col_cc[neighbors[col_cf[j_0]].second] = i_0;
983 }
984 }
985 }
986 }
987
988
989
997 template <int max_n_vertices, typename FU>
998 void
1000 const unsigned int face_dimensionality,
1001 const std::vector<std::shared_ptr<CellTypeBase>> &cell_types,
1002 const std::vector<ReferenceCell> &cell_types_index,
1003 const CRS<unsigned int> &crs,
1004 CRS<unsigned int> &crs_d, // result
1005 CRS<unsigned int> &crs_0, // result
1006 TriaObjectsOrientations &orientations, // result
1007 const FU &second_key_function)
1008 {
1009 const bool compatibility_mode = true;
1010
1011 const std::vector<std::size_t> &cell_ptr = crs.ptr;
1012 const std::vector<unsigned int> &cell_vertices = crs.col;
1013 std::vector<std::size_t> &ptr_d = crs_d.ptr;
1014 std::vector<unsigned int> &col_d = crs_d.col;
1015
1016 // note: we do not pre-allocate memory for these arrays because it turned
1017 // out that counting unique entities is more expensive than push_back().
1018 std::vector<std::size_t> &ptr_0 = crs_0.ptr;
1019 std::vector<unsigned int> &col_0 = crs_0.col;
1020
1021 // clear
1022 ptr_0 = {};
1023 col_0 = {};
1024
1025 unsigned int n_entities = 0;
1026
1027 for (const auto &c : cell_types_index)
1028 n_entities +=
1029 cell_types[static_cast<types::geometric_entity_type>(c)]->n_entities(
1030 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 =
1059 cell_types[static_cast<types::geometric_entity_type>(
1060 cell_types_index[c])];
1061 ptr_d[c + 1] = ptr_d[c] + cell_type->n_entities(face_dimensionality);
1062
1063 // ... collect vertices of cell
1064 const ::ArrayView<const unsigned int> cell_vertice(
1065 cell_vertices.data() + cell_ptr[c], cell_ptr[c + 1] - cell_ptr[c]);
1066
1067 // ... loop over all its entities
1068 for (unsigned int e = 0;
1069 e < cell_type->n_entities(face_dimensionality);
1070 ++e)
1071 {
1072 // ... determine global entity vertices
1073 const auto &local_entity_vertices =
1074 cell_type->vertices_of_entity(face_dimensionality, e);
1075
1076 std::array<unsigned int, max_n_vertices> entity_vertices;
1077 std::fill(entity_vertices.begin(), entity_vertices.end(), 0);
1078
1079 for (unsigned int i = 0; i < local_entity_vertices.size(); ++i)
1080 entity_vertices[i] =
1081 cell_vertice[local_entity_vertices[i]] + offset;
1082
1083 // ... create key
1084 std::array<unsigned int, max_n_vertices> key = entity_vertices;
1085 std::sort(key.begin(), key.end());
1086 keys.emplace_back(key, counter++);
1087
1088 ad_entity_vertices.emplace_back(entity_vertices);
1089
1090 ad_entity_types.emplace_back(
1091 cell_type->type_of_entity(face_dimensionality, e));
1092
1093 if (compatibility_mode)
1094 ad_compatibility.emplace_back(
1095 second_key_function(entity_vertices, cell_type, c, e));
1096 }
1097 }
1098
1099 col_d.resize(keys.size());
1100 orientations.reinit(keys.size());
1101
1102 // step 2: sort according to key so that entities with same key can be
1103 // merged
1104 std::sort(keys.begin(), keys.end());
1105
1106
1107 if (compatibility_mode)
1108 {
1109 unsigned int n_unique_entities = 0;
1110 unsigned int n_unique_entity_vertices = 0;
1111
1112 std::array<unsigned int, max_n_vertices> ref_key, new_key;
1113 std::fill(ref_key.begin(), ref_key.end(), 0);
1114 for (unsigned int i = 0; i < keys.size(); ++i)
1115 {
1116 const auto offset_i = std::get<1>(keys[i]);
1117
1118 if (ref_key != std::get<0>(keys[i]))
1119 {
1120 ref_key = std::get<0>(keys[i]);
1121
1122 ++n_unique_entities;
1123 n_unique_entity_vertices +=
1124 cell_types[static_cast<types::geometric_entity_type>(
1125 ad_entity_types[offset_i])]
1126 ->n_entities(0);
1127
1128 new_key = ad_compatibility[offset_i];
1129 }
1130
1131 std::get<0>(keys[i]) = new_key;
1132 }
1133
1134 std::sort(keys.begin(), keys.end());
1135
1136 ptr_0.reserve(n_unique_entities + 1);
1137 col_0.reserve(n_unique_entity_vertices);
1138 }
1139
1140
1141 std::array<unsigned int, max_n_vertices> ref_key;
1142 std::array<unsigned int, max_n_vertices> ref_indices;
1143 std::fill(ref_key.begin(), ref_key.end(), 0);
1144
1145 unsigned int counter = ::numbers::invalid_unsigned_int;
1146 for (unsigned int i = 0; i < keys.size(); i++)
1147 {
1148 const auto offset_i = std::get<1>(keys[i]);
1149
1150 if (ref_key != std::get<0>(keys[i]))
1151 {
1152 // new key: default orientation is correct
1153 ++counter;
1154 ref_key = std::get<0>(keys[i]);
1155 ref_indices = ad_entity_vertices[offset_i];
1156
1157 ptr_0.push_back(col_0.size());
1158 for (const auto j : ad_entity_vertices[offset_i])
1159 if (j != 0)
1160 col_0.push_back(j - offset);
1161 }
1162 else
1163 {
1164 // previously seen key: set orientation relative to the first
1165 // occurrence
1166 orientations.set_combined_orientation(
1167 offset_i,
1168 ad_entity_types[offset_i]
1169 .template get_combined_orientation<unsigned int>(
1170 make_array_view(ad_entity_vertices[offset_i].begin(),
1171 ad_entity_vertices[offset_i].begin() +
1172 ad_entity_types[offset_i].n_vertices()),
1173 make_array_view(ref_indices.begin(),
1174 ref_indices.begin() +
1175 ad_entity_types[offset_i].n_vertices())));
1176 }
1177 col_d[offset_i] = counter;
1178 }
1179 ptr_0.push_back(col_0.size());
1180 }
1181
1182
1183
1188 template <typename FU>
1189 void
1191 const unsigned int face_dimensionality,
1192 const std::vector<std::shared_ptr<CellTypeBase>> &cell_types,
1193 const std::vector<ReferenceCell> &cell_types_index,
1194 const CRS<unsigned int> &crs,
1195 CRS<unsigned int> &crs_d,
1196 CRS<unsigned int> &crs_0,
1197 TriaObjectsOrientations &orientations,
1198 const FU &second_key_function)
1199 {
1200 std::size_t max_n_vertices = 0;
1201
1202 for (const auto &c : cell_types_index)
1203 {
1204 const auto &cell_type =
1205 cell_types[static_cast<types::geometric_entity_type>(c)];
1206 for (unsigned int e = 0;
1207 e < cell_type->n_entities(face_dimensionality);
1208 ++e)
1209 max_n_vertices = std::max(
1210 max_n_vertices,
1211 cell_type->vertices_of_entity(face_dimensionality, e).size());
1212 }
1213
1214 if (max_n_vertices == 2)
1215 build_face_entities_templated<2>(face_dimensionality,
1216 cell_types,
1217 cell_types_index,
1218 crs,
1219 crs_d,
1220 crs_0,
1221 orientations,
1222 second_key_function);
1223 else if (max_n_vertices == 3)
1224 build_face_entities_templated<3>(face_dimensionality,
1225 cell_types,
1226 cell_types_index,
1227 crs,
1228 crs_d,
1229 crs_0,
1230 orientations,
1231 second_key_function);
1232 else if (max_n_vertices == 4)
1233 build_face_entities_templated<4>(face_dimensionality,
1234 cell_types,
1235 cell_types_index,
1236 crs,
1237 crs_d,
1238 crs_0,
1239 orientations,
1240 second_key_function);
1241 else
1243 }
1244
1245
1246
1254 inline void
1256 const std::vector<std::shared_ptr<CellTypeBase>> &cell_types,
1257 const std::vector<ReferenceCell> &cell_types_index,
1258 const CRS<unsigned int> &con_cv,
1259 const CRS<unsigned int> &con_cl,
1260 const CRS<unsigned int> &con_lv,
1261 const CRS<unsigned int> &con_cq,
1262 const CRS<unsigned int> &con_qv,
1263 const TriaObjectsOrientations &ori_cq,
1264 CRS<unsigned int> &con_ql, // result
1265 TriaObjectsOrientations &ori_ql, // result
1266 std::vector<ReferenceCell> &quad_t_id // result
1267 )
1268 {
1269 // reset output
1270 con_ql.ptr = {};
1271 con_ql.col = {};
1272
1273 con_ql.ptr.resize(con_qv.ptr.size());
1274 con_ql.ptr[0] = 0;
1275
1276 quad_t_id.resize(con_qv.ptr.size() - 1);
1277
1278 // count the number of lines of each face
1279 for (unsigned int c = 0; c < con_cq.ptr.size() - 1; ++c)
1280 {
1281 const auto &cell_type =
1282 cell_types[static_cast<types::geometric_entity_type>(
1283 cell_types_index[c])];
1284
1285 // loop over faces
1286 for (unsigned int f_ = con_cq.ptr[c], f_index = 0;
1287 f_ < con_cq.ptr[c + 1];
1288 ++f_, ++f_index)
1289 {
1290 const unsigned int f = con_cq.col[f_];
1291
1292 con_ql.ptr[f + 1] = cell_type->n_lines_of_surface(f_index);
1293 }
1294 }
1295
1296 // use the counts to determine the offsets -> prefix sum
1297 for (unsigned int i = 0; i < con_ql.ptr.size() - 1; ++i)
1298 con_ql.ptr[i + 1] += con_ql.ptr[i];
1299
1300 // allocate memory
1301 con_ql.col.resize(con_ql.ptr.back());
1302 ori_ql.reinit(con_ql.ptr.back());
1303
1304 // loop over cells
1305 for (unsigned int c = 0; c < con_cq.ptr.size() - 1; ++c)
1306 {
1307 const auto &cell_type =
1308 cell_types[static_cast<types::geometric_entity_type>(
1309 cell_types_index[c])];
1310
1311 // loop over faces
1312 for (unsigned int f_ = con_cq.ptr[c], f_index = 0;
1313 f_ < con_cq.ptr[c + 1];
1314 ++f_, ++f_index)
1315 {
1316 const unsigned int f = con_cq.col[f_];
1317
1318 // only faces with default orientation have to do something
1319 if (ori_cq.get_combined_orientation(f_) !=
1321 continue;
1322
1323 // determine entity type of face
1324 quad_t_id[f] = cell_type->type_of_entity(2, f_index);
1325
1326 // loop over lines
1327 for (unsigned int l = 0;
1328 l < cell_type->n_lines_of_surface(f_index);
1329 ++l)
1330 {
1331 // determine global index of line
1332 const unsigned int local_line_index =
1333 cell_type->nth_line_of_surface(l, f_index);
1334 const unsigned int global_line_index =
1335 con_cl.col[con_cl.ptr[c] + local_line_index];
1336 con_ql.col[con_ql.ptr[f] + l] = global_line_index;
1337
1338 // determine orientation of line
1339 const auto line_vertices_1_ref =
1340 cell_type->vertices_of_nth_line_of_surface(l, f_index);
1341
1342 bool same = true;
1343 for (unsigned int v = 0; v < line_vertices_1_ref.size(); ++v)
1344 if (con_cv.col[con_cv.ptr[c] + line_vertices_1_ref[v]] !=
1345 con_lv.col[con_lv.ptr[global_line_index] + v])
1346 {
1347 same = false;
1348 break;
1349 }
1350
1351 // ... comparison gives orientation
1353 con_ql.ptr[f] + l,
1356 }
1357 }
1358 }
1359 }
1360
1361
1362
1371 template <typename T>
1372 Connectivity<T>
1373 build_connectivity(const unsigned int dim,
1374 const std::vector<std::shared_ptr<CellTypeBase>> &cell_t,
1375 const std::vector<ReferenceCell> &cell_t_id,
1376 const CRS<T> &con_cv)
1377 {
1378 Connectivity<T> connectivity(dim, cell_t_id);
1379
1380 CRS<T> temp1; // needed for 3d
1381
1382 if (dim == 1)
1383 connectivity.entity_to_entities(1, 0) = con_cv;
1384
1385 if (dim == 2 || dim == 3) // build lines
1386 {
1388
1390 1,
1391 cell_t,
1392 connectivity.entity_types(dim),
1393 con_cv,
1394 dim == 2 ? connectivity.entity_to_entities(2, 1) : temp1,
1395 connectivity.entity_to_entities(1, 0),
1396 dim == 2 ? connectivity.entity_orientations(1) : dummy,
1397 [](auto key, const auto &, const auto &, const auto &) {
1398 // to ensure same enumeration as in deal.II
1399 return key;
1400 });
1401 }
1402
1403 if (dim == 3) // build quads
1404 {
1406 2,
1407 cell_t,
1408 connectivity.entity_types(3),
1409 con_cv,
1410 connectivity.entity_to_entities(3, 2),
1411 connectivity.entity_to_entities(2, 0),
1412 connectivity.entity_orientations(2),
1413 [&](auto key, const auto &cell_type, const auto &c, const auto &f) {
1414 // to ensure same enumeration as in deal.II
1415 AssertIndexRange(cell_type->n_lines_of_surface(f),
1416 key.size() + 1);
1417
1418 unsigned int l = 0;
1419
1420 for (; l < cell_type->n_lines_of_surface(f); ++l)
1421 key[l] =
1422 temp1
1423 .col[temp1.ptr[c] + cell_type->nth_line_of_surface(l, f)] +
1424 1 /*offset!*/;
1425
1426 for (; l < key.size(); ++l)
1427 key[l] = 0;
1428
1429 return key;
1430 });
1431
1432 // create connectivity: quad -> line
1433 build_intersection(cell_t,
1434 connectivity.entity_types(3),
1435 con_cv,
1436 temp1,
1437 connectivity.entity_to_entities(1, 0),
1438 connectivity.entity_to_entities(3, 2),
1439 connectivity.entity_to_entities(2, 0),
1440 connectivity.entity_orientations(2),
1441 connectivity.entity_to_entities(2, 1),
1442 connectivity.entity_orientations(1),
1443 connectivity.entity_types(2));
1444 }
1445
1446 // determine neighbors
1447 determine_neighbors(connectivity.entity_to_entities(dim, dim - 1),
1448 connectivity.entity_to_entities(dim, dim));
1449
1450 return connectivity;
1451 }
1452
1453
1454
1458 template <typename T, int dim>
1459 Connectivity<T>
1460 build_connectivity(const std::vector<CellData<dim>> &cells)
1461 {
1462 AssertThrow(cells.size() > 0, ExcMessage("No cells have been provided!"));
1463
1464 // vector of possible cell entity types
1465 std::vector<std::shared_ptr<CellTypeBase>> cell_types_impl(8);
1466
1467 cell_types_impl[static_cast<types::geometric_entity_type>(
1468 ReferenceCells::Line)] = std::make_shared<CellTypeLine>();
1469 cell_types_impl[static_cast<types::geometric_entity_type>(
1470 ReferenceCells::Triangle)] = std::make_shared<CellTypeTriangle>();
1471 cell_types_impl[static_cast<types::geometric_entity_type>(
1473 std::make_shared<CellTypeQuadrilateral>();
1474 cell_types_impl[static_cast<types::geometric_entity_type>(
1475 ReferenceCells::Tetrahedron)] = std::make_shared<CellTypeTetrahedron>();
1476 cell_types_impl[static_cast<types::geometric_entity_type>(
1477 ReferenceCells::Pyramid)] = std::make_shared<CellTypePyramid>();
1478 cell_types_impl[static_cast<types::geometric_entity_type>(
1479 ReferenceCells::Wedge)] = std::make_shared<CellTypeWedge>();
1480 cell_types_impl[static_cast<types::geometric_entity_type>(
1481 ReferenceCells::Hexahedron)] = std::make_shared<CellTypeHexahedron>();
1482
1483 // determine cell types and process vertices
1484 std::vector<T> cell_vertices;
1485 cell_vertices.reserve(
1486 std::accumulate(cells.begin(),
1487 cells.end(),
1488 0,
1489 [](const auto &result, const auto &cell) {
1490 return result + cell.vertices.size();
1491 }));
1492
1493 std::vector<std::size_t> cell_vertices_ptr;
1494 cell_vertices_ptr.reserve(cells.size() + 1);
1495 cell_vertices_ptr.push_back(0);
1496
1497 std::vector<ReferenceCell> cell_types_indices;
1498 cell_types_indices.reserve(cells.size());
1499
1500 // loop over cells and create CRS
1501 for (const auto &cell : cells)
1502 {
1503#ifdef DEBUG
1504 auto vertices_unique = cell.vertices;
1505 std::sort(vertices_unique.begin(), vertices_unique.end());
1506 vertices_unique.erase(std::unique(vertices_unique.begin(),
1507 vertices_unique.end()),
1508 vertices_unique.end());
1509
1510 Assert(vertices_unique.size() == cell.vertices.size(),
1511 ExcMessage(
1512 "The definition of a cell refers to the same vertex several "
1513 "times. This is not possible. A common reason is that "
1514 "CellData::vertices has a size that does not match the "
1515 "size expected from the reference cell. Please resize "
1516 "CellData::vertices or use the appropriate constructor of "
1517 "CellData."));
1518#endif
1519
1520 const ReferenceCell reference_cell =
1521 ReferenceCell::n_vertices_to_type(dim, cell.vertices.size());
1522
1523 Assert(reference_cell != ReferenceCells::Invalid,
1526 reference_cell),
1527 cell_types_impl.size());
1528 Assert(cell_types_impl[static_cast<types::geometric_entity_type>(
1529 reference_cell)]
1530 .get() != nullptr,
1532
1533 cell_types_indices.push_back(reference_cell);
1534
1535 // create CRS of vertices (to remove template argument dim)
1536 for (const auto &vertex : cell.vertices)
1537 cell_vertices.push_back(vertex);
1538
1539 cell_vertices_ptr.push_back(cell_vertices.size());
1540 }
1541
1542 // do the actual work
1543 return build_connectivity<T>(dim,
1544 cell_types_impl,
1545 cell_types_indices,
1546 {cell_vertices_ptr, cell_vertices});
1547 }
1548 } // namespace TriangulationImplementation
1549} // namespace internal
1550
1551
1553
1554#endif
ArrayView< typename std::remove_reference< typename std::iterator_traits< Iterator >::reference >::type, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition array_view.h:949
static constexpr unsigned char default_combined_face_orientation()
static constexpr unsigned char reversed_combined_line_orientation()
static ReferenceCell n_vertices_to_type(const int dim, const unsigned int n_vertices)
void set_combined_orientation(const unsigned int object, const unsigned char value)
unsigned char get_combined_orientation(const unsigned int object) const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
Point< 2 > first
Definition grid_out.cc:4623
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)
static const unsigned int invalid_unsigned_int
Definition types.h:220
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
std::uint8_t geometric_entity_type
Definition types.h:172
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