Loading [MathJax]/extensions/TeX/newcommand.js
 Reference documentation for deal.II version 9.3.3
\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\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
connectivity.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2020 - 2021 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
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 {
52 Assert(false, ExcNotImplemented());
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 {
64 Assert(false, ExcNotImplemented());
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 {
77 Assert(false, ExcNotImplemented());
78 (void)d;
79 (void)e;
80
82 }
83
87 virtual unsigned int
88 n_lines_of_surface(const unsigned int face) const
89 {
90 Assert(false, ExcNotImplemented());
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 {
103 Assert(false, ExcNotImplemented());
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 {
117 Assert(false, ExcNotImplemented());
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
149 Assert(false, ExcNotImplemented());
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
162 Assert(false, ExcNotImplemented());
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
180 struct CellTypeTri : public CellTypeBase
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
203 Assert(false, ExcNotImplemented());
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
219 Assert(false, ExcNotImplemented());
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
260 Assert(false, ExcNotImplemented());
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
276 Assert(false, ExcNotImplemented());
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
294 struct CellTypeTet : public CellTypeBase
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
325 Assert(false, ExcNotImplemented());
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
344 Assert(false, ExcNotImplemented());
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
436 Assert(false, ExcNotImplemented());
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
457 Assert(false, ExcNotImplemented());
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
562 Assert(false, ExcNotImplemented());
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
583 Assert(false, ExcNotImplemented());
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
642 struct CellTypeHex : public CellTypeBase
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
690 Assert(false, ExcNotImplemented());
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
709 Assert(false, ExcNotImplemented());
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
830 inline std::vector<unsigned char> &
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 std::vector<unsigned char> &
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
892 Assert(false, ExcNotImplemented());
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
911 Assert(false, ExcNotImplemented());
912
913 return cell_entities;
914 }
915
916 private:
917 const unsigned int dim;
918 std::vector<::ReferenceCell> cell_types;
919
921
922 std::vector<unsigned char> line_orientation;
923
926
927 std::vector<unsigned char> quad_orientation;
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 key_length, typename FU>
998 void
1000 const unsigned int d,
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 std::vector<unsigned char> & 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 d);
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, key_length>, unsigned int>>
1039 keys; // key (sorted vertices), cell-entity index
1040
1041 std::vector<std::array<unsigned int, key_length>> ad_entity_vertices;
1042 std::vector<::ReferenceCell> ad_entity_types;
1043 std::vector<std::array<unsigned int, key_length>> 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(d);
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; e < cell_type->n_entities(d); e++)
1069 {
1070 // ... determine global entity vertices
1071 const auto &local_entity_vertices =
1072 cell_type->vertices_of_entity(d, e);
1073
1074 std::array<unsigned int, key_length> 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, key_length> 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(cell_type->type_of_entity(d, e));
1089
1090 if (compatibility_mode)
1091 ad_compatibility.emplace_back(
1092 second_key_function(entity_vertices, cell_type, c, e));
1093 }
1094 }
1095
1096 col_d.resize(keys.size());
1097 orientations.resize(keys.size());
1098
1099 // step 2: sort according to key so that entities with same key can be
1100 // merged
1101 std::sort(keys.begin(), keys.end());
1102
1103
1104 if (compatibility_mode)
1105 {
1106 unsigned int n_unique_entities = 0;
1107 unsigned int n_unique_entity_vertices = 0;
1108
1109 std::array<unsigned int, key_length> ref_key, new_key;
1110 std::fill(ref_key.begin(), ref_key.end(), 0);
1111 for (unsigned int i = 0; i < keys.size(); ++i)
1112 {
1113 const auto offset_i = std::get<1>(keys[i]);
1114
1115 if (ref_key != std::get<0>(keys[i]))
1116 {
1117 ref_key = std::get<0>(keys[i]);
1118
1119 n_unique_entities++;
1120 n_unique_entity_vertices +=
1121 cell_types[static_cast<types::geometric_entity_type>(
1122 ad_entity_types[offset_i])]
1123 ->n_entities(0);
1124
1125 new_key = ad_compatibility[offset_i];
1126 }
1127
1128 std::get<0>(keys[i]) = new_key;
1129 }
1130
1131 std::sort(keys.begin(), keys.end());
1132
1133 ptr_0.reserve(n_unique_entities);
1134 col_0.reserve(n_unique_entity_vertices);
1135 }
1136
1137
1138 std::array<unsigned int, key_length> ref_key;
1139 std::array<unsigned int, key_length> ref_indices;
1140 std::fill(ref_key.begin(), ref_key.end(), 0);
1141
1142 for (unsigned int i = 0, counter = ::numbers::invalid_unsigned_int;
1143 i < keys.size();
1144 i++)
1145 {
1146 const auto offset_i = std::get<1>(keys[i]);
1147
1148 if (ref_key != std::get<0>(keys[i]))
1149 {
1150 // new key
1151 counter++;
1152 ref_key = std::get<0>(keys[i]);
1153 ref_indices = ad_entity_vertices[offset_i];
1154
1155 ptr_0.push_back(col_0.size());
1156 for (const auto j : ad_entity_vertices[offset_i])
1157 if (j != 0)
1158 col_0.push_back(j - offset);
1159
1160 // take its orientation as default
1161 col_d[offset_i] = counter;
1162 orientations[offset_i] = 1;
1163 }
1164 else
1165 {
1166 col_d[offset_i] = counter;
1167 orientations[offset_i] =
1168 ad_entity_types[offset_i].compute_orientation(
1169 ad_entity_vertices[offset_i], ref_indices);
1170 }
1171 }
1172 ptr_0.push_back(col_0.size());
1173 }
1174
1175
1176
1181 template <typename FU>
1182 void
1183 build_entity(const unsigned int d,
1184 const std::vector<std::shared_ptr<CellTypeBase>> &cell_types,
1185 const std::vector<::ReferenceCell> &cell_types_index,
1186 const CRS<unsigned int> & crs,
1187 CRS<unsigned int> & crs_d,
1188 CRS<unsigned int> & crs_0,
1189 std::vector<unsigned char> & orientations,
1190 const FU & second_key_function)
1191 {
1192 std::size_t key_length = 0;
1193
1194 for (const auto &c : cell_types_index)
1195 {
1196 const auto &cell_type =
1197 cell_types[static_cast<types::geometric_entity_type>(c)];
1198 for (unsigned int e = 0; e < cell_type->n_entities(d); e++)
1199 key_length =
1200 std::max(key_length, cell_type->vertices_of_entity(d, e).size());
1201 }
1202
1203 if (key_length == 2)
1204 build_entity_templated<2>(d,
1205 cell_types,
1206 cell_types_index,
1207 crs,
1208 crs_d,
1209 crs_0,
1210 orientations,
1211 second_key_function);
1212 else if (key_length == 3)
1213 build_entity_templated<3>(d,
1214 cell_types,
1215 cell_types_index,
1216 crs,
1217 crs_d,
1218 crs_0,
1219 orientations,
1220 second_key_function);
1221 else if (key_length == 4)
1222 build_entity_templated<4>(d,
1223 cell_types,
1224 cell_types_index,
1225 crs,
1226 crs_d,
1227 crs_0,
1228 orientations,
1229 second_key_function);
1230 else
1232 }
1233
1234
1235
1243 inline void
1245 const std::vector<std::shared_ptr<CellTypeBase>> &cell_types,
1246 const std::vector<::ReferenceCell> & cell_types_index,
1247 const CRS<unsigned int> & con_cv,
1248 const CRS<unsigned int> & con_cl,
1249 const CRS<unsigned int> & con_lv,
1250 const CRS<unsigned int> & con_cq,
1251 const CRS<unsigned int> & con_qv,
1252 const std::vector<unsigned char> & ori_cq,
1253 CRS<unsigned int> & con_ql, // result
1254 std::vector<unsigned char> & ori_ql, // result
1255 std::vector<::ReferenceCell> & quad_t_id // result
1256 )
1257 {
1258 // reset output
1259 ori_ql = {};
1260 con_ql.ptr = {};
1261 con_ql.col = {};
1262
1263 con_ql.ptr.resize(con_qv.ptr.size());
1264 con_ql.ptr[0] = 0;
1265
1266 quad_t_id.resize(con_qv.ptr.size() - 1);
1267
1268 // count the number of lines of each face
1269 for (unsigned int c = 0; c < con_cq.ptr.size() - 1; ++c)
1270 {
1271 const auto &cell_type =
1272 cell_types[static_cast<types::geometric_entity_type>(
1273 cell_types_index[c])];
1274
1275 // loop over faces
1276 for (unsigned int f_ = con_cq.ptr[c], f_index = 0;
1277 f_ < con_cq.ptr[c + 1];
1278 ++f_, ++f_index)
1279 {
1280 const unsigned int f = con_cq.col[f_];
1281
1282 con_ql.ptr[f + 1] = cell_type->n_lines_of_surface(f_index);
1283 }
1284 }
1285
1286 // use the counts to determine the offsets -> prefix sum
1287 for (unsigned int i = 0; i < con_ql.ptr.size() - 1; ++i)
1288 con_ql.ptr[i + 1] += con_ql.ptr[i];
1289
1290 // allocate memory
1291 con_ql.col.resize(con_ql.ptr.back());
1292 ori_ql.resize(con_ql.ptr.back());
1293
1294 // loop over cells
1295 for (unsigned int c = 0; c < con_cq.ptr.size() - 1; ++c)
1296 {
1297 const auto &cell_type =
1298 cell_types[static_cast<types::geometric_entity_type>(
1299 cell_types_index[c])];
1300
1301 // loop over faces
1302 for (unsigned int f_ = con_cq.ptr[c], f_index = 0;
1303 f_ < con_cq.ptr[c + 1];
1304 ++f_, ++f_index)
1305 {
1306 const unsigned int f = con_cq.col[f_];
1307
1308 // only faces with default orientation have to do something
1309 if (ori_cq[f_] != 1)
1310 continue;
1311
1312 // determine entity type of face
1313 quad_t_id[f] = cell_type->type_of_entity(2, f_index);
1314
1315 // loop over lines
1316 for (unsigned int l = 0;
1317 l < cell_type->n_lines_of_surface(f_index);
1318 ++l)
1319 {
1320 // determine global index of line
1321 const unsigned int local_line_index =
1322 cell_type->nth_line_of_surface(l, f_index);
1323 const unsigned int global_line_index =
1324 con_cl.col[con_cl.ptr[c] + local_line_index];
1325 con_ql.col[con_ql.ptr[f] + l] = global_line_index;
1326
1327 // determine orientation of line
1328 const auto line_vertices_1_ref =
1329 cell_type->vertices_of_nth_line_of_surface(l, f_index);
1330
1331 bool same = true;
1332 for (unsigned int v = 0; v < line_vertices_1_ref.size(); ++v)
1333 if (con_cv.col[con_cv.ptr[c] + line_vertices_1_ref[v]] !=
1334 con_lv.col[con_lv.ptr[global_line_index] + v])
1335 {
1336 same = false;
1337 break;
1338 }
1339
1340 // ... comparison gives orientation
1341 ori_ql[con_ql.ptr[f] + l] = (same ? 1 : 0);
1342 }
1343 }
1344 }
1345 }
1346
1347
1348
1357 template <typename T>
1358 Connectivity<T>
1359 build_connectivity(const unsigned int dim,
1360 const std::vector<std::shared_ptr<CellTypeBase>> &cell_t,
1361 const std::vector<::ReferenceCell> &cell_t_id,
1362 const CRS<T> & con_cv)
1363 {
1364 Connectivity<T> connectivity(dim, cell_t_id);
1365
1366 CRS<T> temp1; // needed for 3D
1367
1368 if (dim == 1)
1369 connectivity.entity_to_entities(1, 0) = con_cv;
1370
1371 if (dim == 2 || dim == 3) // build lines
1372 {
1373 std::vector<unsigned char> dummy;
1374
1375 build_entity(1,
1376 cell_t,
1377 connectivity.entity_types(dim),
1378 con_cv,
1379 dim == 2 ? connectivity.entity_to_entities(2, 1) : temp1,
1380 connectivity.entity_to_entities(1, 0),
1381 dim == 2 ? connectivity.entity_orientations(1) : dummy,
1382 [](auto key, const auto &, const auto &, const auto &) {
1383 // to ensure same enumeration as in deal.II
1384 return key;
1385 });
1386 }
1387
1388 if (dim == 3) // build quads
1389 {
1391 2,
1392 cell_t,
1393 connectivity.entity_types(3),
1394 con_cv,
1395 connectivity.entity_to_entities(3, 2),
1396 connectivity.entity_to_entities(2, 0),
1397 connectivity.entity_orientations(2),
1398 [&](auto key, const auto &cell_type, const auto &c, const auto &f) {
1399 // to ensure same enumeration as in deal.II
1400 AssertIndexRange(cell_type->n_lines_of_surface(f),
1401 key.size() + 1);
1402
1403 unsigned int l = 0;
1404
1405 for (; l < cell_type->n_lines_of_surface(f); ++l)
1406 key[l] =
1407 temp1
1408 .col[temp1.ptr[c] + cell_type->nth_line_of_surface(l, f)] +
1409 1 /*offset!*/;
1410
1411 for (; l < key.size(); ++l)
1412 key[l] = 0;
1413
1414 return key;
1415 });
1416
1417 // create connectivity: quad -> line
1418 build_intersection(cell_t,
1419 connectivity.entity_types(3),
1420 con_cv,
1421 temp1,
1422 connectivity.entity_to_entities(1, 0),
1423 connectivity.entity_to_entities(3, 2),
1424 connectivity.entity_to_entities(2, 0),
1425 connectivity.entity_orientations(2),
1426 connectivity.entity_to_entities(2, 1),
1427 connectivity.entity_orientations(1),
1428 connectivity.entity_types(2));
1429 }
1430
1431 // determine neighbors
1432 determine_neighbors(connectivity.entity_to_entities(dim, dim - 1),
1433 connectivity.entity_to_entities(dim, dim));
1434
1435 return connectivity;
1436 }
1437
1438
1439
1443 template <typename T, int dim>
1444 Connectivity<T>
1445 build_connectivity(const std::vector<CellData<dim>> &cells)
1446 {
1447 AssertThrow(cells.size() > 0, ExcMessage("No cells have been provided!"));
1448
1449 // vector of possible cell entity types
1450 std::vector<std::shared_ptr<CellTypeBase>> cell_types_impl(8);
1451
1452 cell_types_impl[static_cast<types::geometric_entity_type>(
1454 .reset(new CellTypeLine());
1455 cell_types_impl[static_cast<types::geometric_entity_type>(
1457 .reset(new CellTypeTri());
1458 cell_types_impl[static_cast<types::geometric_entity_type>(
1460 .reset(new CellTypeQuad());
1461 cell_types_impl[static_cast<types::geometric_entity_type>(
1463 .reset(new CellTypeTet());
1464 cell_types_impl[static_cast<types::geometric_entity_type>(
1466 .reset(new CellTypePyramid());
1467 cell_types_impl[static_cast<types::geometric_entity_type>(
1469 .reset(new CellTypeWedge());
1470 cell_types_impl[static_cast<types::geometric_entity_type>(
1472 .reset(new CellTypeHex());
1473
1474 // determine cell types and process vertices
1475 std::vector<T> cell_vertices;
1476 cell_vertices.reserve(
1477 std::accumulate(cells.begin(),
1478 cells.end(),
1479 0,
1480 [](const auto &result, const auto &cell) {
1481 return result + cell.vertices.size();
1482 }));
1483
1484 std::vector<std::size_t> cell_vertices_ptr;
1485 cell_vertices_ptr.reserve(cells.size() + 1);
1486 cell_vertices_ptr.push_back(0);
1487
1488 std::vector<::ReferenceCell> cell_types_indices;
1489 cell_types_indices.reserve(cells.size());
1490
1491 // loop over cells and create CRS
1492 for (const auto &cell : cells)
1493 {
1494 const ::ReferenceCell reference_cell =
1496 cell.vertices.size());
1497
1502 cell_types_impl.size());
1503 Assert(cell_types_impl[static_cast<types::geometric_entity_type>(
1505 .get() != nullptr,
1507
1508 cell_types_indices.push_back(reference_cell);
1509
1510 // create CRS of vertices (to remove template argument dim)
1511 for (const auto &vertex : cell.vertices)
1512 cell_vertices.push_back(vertex);
1513
1514 cell_vertices_ptr.push_back(cell_vertices.size());
1515 }
1516
1517 // do the actual work
1518 return build_connectivity<T>(dim,
1519 cell_types_impl,
1520 cell_types_indices,
1521 {cell_vertices_ptr, cell_vertices});
1522 }
1523 } // namespace TriangulationImplementation
1524} // namespace internal
1525
1526
1528
1529#endif
static ReferenceCell n_vertices_to_type(const int dim, const unsigned int n_vertices)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:402
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
Point< 2 > first
Definition: grid_out.cc:4587
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
Definition: exceptions.h:1465
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1575
graph_traits< Graph >::vertex_descriptor Vertex
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
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 Line
void determine_neighbors(const CRS< T > &con_cf, CRS< T > &con_cc)
Definition: connectivity.h:946
void build_entity_templated(const unsigned int d, 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, std::vector< unsigned char > &orientations, const FU &second_key_function)
Definition: connectivity.h:999
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 std::vector< unsigned char > &ori_cq, CRS< unsigned int > &con_ql, std::vector< unsigned char > &ori_ql, std::vector<::ReferenceCell > &quad_t_id)
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_entity(const unsigned int d, 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, std::vector< unsigned char > &orientations, const FU &second_key_function)
static const unsigned int invalid_unsigned_int
Definition: types.h:196
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
std::uint8_t geometric_entity_type
Definition: types.h:157
CRS(const std::vector< std::size_t > &ptr, const std::vector< T > &col)
Definition: connectivity.h:779
virtual unsigned int n_entities(const unsigned int d) const
Definition: connectivity.h:50
virtual ::ArrayView< const unsigned int > vertices_of_entity(const unsigned int d, const unsigned int e) const
Definition: connectivity.h:62
virtual ::ReferenceCell type_of_entity(const unsigned int d, const unsigned int e) const
Definition: connectivity.h:75
virtual const std::array< unsigned int, 2 > & vertices_of_nth_line_of_surface(const unsigned int line, const unsigned int face) const
Definition: connectivity.h:114
virtual unsigned int n_lines_of_surface(const unsigned int face) const
Definition: connectivity.h:88
virtual unsigned int nth_line_of_surface(const unsigned int line, const unsigned int face) const
Definition: connectivity.h:100
virtual ::ReferenceCell type_of_entity(const unsigned int d, const unsigned int e) const override
Definition: connectivity.h:696
::ArrayView< const unsigned int > vertices_of_entity(const unsigned int d, const unsigned int e) const override
Definition: connectivity.h:645
unsigned int n_entities(const unsigned int d) const override
Definition: connectivity.h:715
const std::array< unsigned int, 2 > & vertices_of_nth_line_of_surface(const unsigned int line, const unsigned int face) const override
Definition: connectivity.h:744
unsigned int nth_line_of_surface(const unsigned int line, const unsigned int face) const override
Definition: connectivity.h:729
unsigned int n_lines_of_surface(const unsigned int surface) const override
Definition: connectivity.h:722
::ReferenceCell type_of_entity(const unsigned int d, const unsigned int e) const override
Definition: connectivity.h:155
::ArrayView< const unsigned int > vertices_of_entity(const unsigned int d, const unsigned int e) const override
Definition: connectivity.h:135
unsigned int n_entities(const unsigned int d) const override
Definition: connectivity.h:168
unsigned int n_lines_of_surface(const unsigned int surface) const override
Definition: connectivity.h:470
virtual ::ReferenceCell type_of_entity(const unsigned int d, const unsigned int e) const override
Definition: connectivity.h:442
const std::array< unsigned int, 2 > & vertices_of_nth_line_of_surface(const unsigned int line, const unsigned int face) const override
Definition: connectivity.h:493
::ArrayView< const unsigned int > vertices_of_entity(const unsigned int d, const unsigned int e) const override
Definition: connectivity.h:395
unsigned int n_entities(const unsigned int d) const override
Definition: connectivity.h:463
unsigned int nth_line_of_surface(const unsigned int line, const unsigned int face) const override
Definition: connectivity.h:479
unsigned int n_entities(const unsigned int d) const override
Definition: connectivity.h:282
::ArrayView< const unsigned int > vertices_of_entity(const unsigned int d, const unsigned int e) const override
Definition: connectivity.h:240
virtual ::ReferenceCell type_of_entity(const unsigned int d, const unsigned int e) const override
Definition: connectivity.h:266
unsigned int nth_line_of_surface(const unsigned int line, const unsigned int face) const override
Definition: connectivity.h:364
unsigned int n_lines_of_surface(const unsigned int line) const override
Definition: connectivity.h:357
const std::array< unsigned int, 2 > & vertices_of_nth_line_of_surface(const unsigned int line, const unsigned int face) const override
Definition: connectivity.h:374
virtual ::ReferenceCell type_of_entity(const unsigned int d, const unsigned int e) const override
Definition: connectivity.h:331
unsigned int n_entities(const unsigned int d) const override
Definition: connectivity.h:350
::ArrayView< const unsigned int > vertices_of_entity(const unsigned int d, const unsigned int e) const override
Definition: connectivity.h:297
unsigned int n_entities(const unsigned int d) const override
Definition: connectivity.h:225
::ArrayView< const unsigned int > vertices_of_entity(const unsigned int d, const unsigned int e) const override
Definition: connectivity.h:183
virtual ::ReferenceCell type_of_entity(const unsigned int d, const unsigned int e) const override
Definition: connectivity.h:209
unsigned int n_entities(const unsigned int d) const override
Definition: connectivity.h:589
unsigned int n_lines_of_surface(const unsigned int surface) const override
Definition: connectivity.h:596
::ArrayView< const unsigned int > vertices_of_entity(const unsigned int d, const unsigned int e) const override
Definition: connectivity.h:517
virtual ::ReferenceCell type_of_entity(const unsigned int d, const unsigned int e) const override
Definition: connectivity.h:568
unsigned int nth_line_of_surface(const unsigned int line, const unsigned int face) const override
Definition: connectivity.h:605
const std::array< unsigned int, 2 > & vertices_of_nth_line_of_surface(const unsigned int line, const unsigned int face) const override
Definition: connectivity.h:621
const CRS< T > & entity_to_entities(const unsigned int from, const unsigned int to) const
Definition: connectivity.h:898
std::vector<::ReferenceCell > & entity_types(const unsigned int structdim)
Definition: connectivity.h:853
Connectivity(const unsigned int dim, const std::vector<::ReferenceCell > &cell_types)
Definition: connectivity.h:824
CRS< T > & entity_to_entities(const unsigned int from, const unsigned int to)
Definition: connectivity.h:879
std::vector< unsigned char > & entity_orientations(const unsigned int structdim)
Definition: connectivity.h:831
const std::vector< unsigned char > & entity_orientations(const unsigned int structdim) const
Definition: connectivity.h:842
const std::vector<::ReferenceCell > & entity_types(const unsigned int structdim) const
Definition: connectivity.h:866