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