deal.II version GIT relicensing-2581-gae2745de1b 2025-02-07 22:30:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
hanging_nodes_internal.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2018 - 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_hanging_nodes_internal_h
16#define dealii_hanging_nodes_internal_h
17
18#include <deal.II/base/config.h>
19
22
24
25#include <deal.II/fe/fe_q.h>
27#include <deal.II/fe/fe_tools.h>
28
30
32
33namespace internal
34{
35 namespace MatrixFreeFunctions
36 {
50 enum class ConstraintKinds : std::uint16_t
51 {
52 // default: unconstrained cell
53 unconstrained = 0,
54
55 // subcell
56 subcell_x = 1 << 0,
57 subcell_y = 1 << 1,
58 subcell_z = 1 << 2,
59
60 // face is constrained
61 face_x = 1 << 3,
62 face_y = 1 << 4,
63 face_z = 1 << 5,
64
65 // edge is constrained
66 edge_x = 1 << 6,
67 edge_y = 1 << 7,
68 edge_z = 1 << 8
69 };
70
71
72
76 using compressed_constraint_kind = std::uint8_t;
77
78
79
85
86
87
91 inline bool
92 check(const ConstraintKinds kind_in, const unsigned int dim)
93 {
94 const std::uint16_t kind = static_cast<std::uint16_t>(kind_in);
95 const std::uint16_t subcell = (kind >> 0) & 7;
96 const std::uint16_t face = (kind >> 3) & 7;
97 const std::uint16_t edge = (kind >> 6) & 7;
98
99 if ((kind >> 9) > 0)
100 return false;
101
102 if (dim == 2)
103 {
104 if (edge > 0)
105 return false; // in 2d there are no edge constraints
106
107 if (subcell == 0 && face == 0)
108 return true; // no constraints
109 else if (0 < face)
110 return true; // at least one face is constrained
111 }
112 else if (dim == 3)
113 {
114 if (subcell == 0 && face == 0 && edge == 0)
115 return true; // no constraints
116 else if (0 < face && edge == 0)
117 return true; // at least one face is constrained
118 else if (0 == face && 0 < edge)
119 return true; // at least one edge is constrained
120 else if ((face == edge) && (face == 1 || face == 2 || face == 4))
121 return true; // one face and its orthogonal edge is constrained
122 }
123
124 return false;
125 }
126
127
128
134 compress(const ConstraintKinds kind_in, const unsigned int dim)
135 {
136 Assert(check(kind_in, dim), ExcInternalError());
137
138 if (dim == 2)
139 return static_cast<compressed_constraint_kind>(kind_in);
140
141 if (kind_in == ConstraintKinds::unconstrained)
143
144 const std::uint16_t kind = static_cast<std::uint16_t>(kind_in);
145 const std::uint16_t subcell = (kind >> 0) & 7;
146 const std::uint16_t face = (kind >> 3) & 7;
147 const std::uint16_t edge = (kind >> 6) & 7;
148
149 return subcell + ((face > 0) << 3) + ((edge > 0) << 4) +
150 (std::max(face, edge) << 5);
151 }
152
153
154
159 inline ConstraintKinds
160 decompress(const compressed_constraint_kind kind_in, const unsigned int dim)
161 {
162 if (dim == 2)
163 return static_cast<ConstraintKinds>(kind_in);
164
167
168 const std::uint16_t subcell = (kind_in >> 0) & 7;
169 const std::uint16_t flag_0 = (kind_in >> 3) & 3;
170 const std::uint16_t flag_1 = (kind_in >> 5) & 7;
171
172 const auto result = static_cast<ConstraintKinds>(
173 subcell + (((flag_0 & 0b01) ? flag_1 : 0) << 3) +
174 (((flag_0 & 0b10) ? flag_1 : 0) << 6));
175
176 Assert(check(result, dim), ExcInternalError());
177
178 return result;
179 }
180
181
182
186 inline std::size_t
188 {
189 return sizeof(ConstraintKinds);
190 }
191
192
193
203 {
204 return static_cast<ConstraintKinds>(static_cast<std::uint16_t>(f1) |
205 static_cast<std::uint16_t>(f2));
206 }
207
208
209
216 {
217 f1 = f1 | f2;
218 return f1;
219 }
220
221
222
226 DEAL_II_HOST_DEVICE inline bool
228 {
229 return static_cast<std::uint16_t>(f1) != static_cast<std::uint16_t>(f2);
230 }
231
232
233
239 operator<(const ConstraintKinds f1, const ConstraintKinds f2)
240 {
241 return static_cast<std::uint16_t>(f1) < static_cast<std::uint16_t>(f2);
242 }
243
244
245
251 {
252 return static_cast<ConstraintKinds>(static_cast<std::uint16_t>(f1) &
253 static_cast<std::uint16_t>(f2));
254 }
255
256
257
264 template <int dim>
266 {
267 public:
271 HangingNodes(const Triangulation<dim> &triangualtion);
272
276 std::size_t
277 memory_consumption() const;
278
282 template <typename CellIterator>
283 bool
285 const CellIterator &cell,
286 const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner,
287 const std::vector<std::vector<unsigned int>> &lexicographic_mapping,
288 std::vector<types::global_dof_index> &dof_indices,
289 const ArrayView<ConstraintKinds> &mask) const;
290
295 static std::vector<std::vector<bool>>
296 compute_supported_components(const ::hp::FECollection<dim> &fe);
297
301 template <typename CellIterator>
303 compute_refinement_configuration(const CellIterator &cell) const;
304
309 template <typename CellIterator>
310 void
312 const CellIterator &cell,
313 const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner,
314 const std::vector<std::vector<unsigned int>> &lexicographic_mapping,
315 const std::vector<std::vector<bool>> &component_mask,
316 const ConstraintKinds &refinement_configuration,
317 std::vector<types::global_dof_index> &dof_indices) const;
318
319 private:
323 void
325
326 void
327 rotate_subface_index(int times, unsigned int &subface_index) const;
328
329 void
330 orient_face(const types::geometric_orientation combined_orientation,
331 const unsigned int n_dofs_1d,
332 std::vector<types::global_dof_index> &dofs) const;
333
334 unsigned int
335 line_dof_idx(int local_line,
336 unsigned int dof,
337 unsigned int n_dofs_1d) const;
338
339 void
340 transpose_subface_index(unsigned int &subface) const;
341
342 std::vector<
343 boost::container::small_vector<std::array<unsigned int, 3>, 6>>
345
346 const ::ndarray<unsigned int, 3, 2, 2> local_lines = {
347 {{{{{7, 3}}, {{6, 2}}}},
348 {{{{5, 1}}, {{4, 0}}}},
349 {{{{11, 9}}, {{10, 8}}}}}};
350 };
351
352
353
354 template <int dim>
357 {
358 // Set up line-to-cell mapping for edge constraints (only if dim = 3 and
359 // for pure hex meshes)
360 if (triangulation.all_reference_cells_are_hyper_cube())
361 setup_line_to_cell(triangulation);
362 }
363
364
365
366 template <int dim>
367 inline std::size_t
369 {
370 std::size_t size = 0;
371 for (const auto &a : line_to_cells)
372 size +=
373 (a.capacity() > 6 ? a.capacity() : 0) * sizeof(a[0]) + sizeof(a);
374 return size;
375 }
376
377
378
379 template <int dim>
380 inline void
386
387
388
389 template <>
390 inline void
392 {
393 // Check if we there are no hanging nodes on the current MPI process,
394 // which we do by checking if the second finest level holds no active
395 // non-artificial cell
396 if (triangulation.n_levels() <= 1 ||
397 std::none_of(triangulation.begin_active(triangulation.n_levels() - 2),
398 triangulation.end_active(triangulation.n_levels() - 2),
399 [](const CellAccessor<3, 3> &cell) {
400 return !cell.is_artificial();
401 }))
402 return;
403
404 const unsigned int n_raw_lines = triangulation.n_raw_lines();
405 this->line_to_cells.resize(n_raw_lines);
406
407 // In 3d, we can have DoFs on only an edge being constrained (e.g. in a
408 // cartesian 2x2x2 grid, where only the upper left 2 cells are refined).
409 // This sets up a helper data structure in the form of a mapping from
410 // edges (i.e. lines) to neighboring cells.
411
412 // Mapping from an edge to which children that share that edge.
413 const unsigned int line_to_children[12][2] = {{0, 2},
414 {1, 3},
415 {0, 1},
416 {2, 3},
417 {4, 6},
418 {5, 7},
419 {4, 5},
420 {6, 7},
421 {0, 4},
422 {1, 5},
423 {2, 6},
424 {3, 7}};
425
426 std::vector<
427 boost::container::small_vector<std::array<unsigned int, 3>, 6>>
428 line_to_inactive_cells(n_raw_lines);
429
430 // First add active and inactive cells to their lines:
431 for (const auto &cell : triangulation.cell_iterators())
432 {
433 const unsigned int cell_level = cell->level();
434 const unsigned int cell_index = cell->index();
435 for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_cell;
436 ++line)
437 {
438 const unsigned int line_idx = cell->line_index(line);
439 if (cell->is_active())
440 line_to_cells[line_idx].push_back(
441 {{cell_level, cell_index, line}});
442 else
443 line_to_inactive_cells[line_idx].push_back(
444 {{cell_level, cell_index, line}});
445 }
446 }
447
448 // Now, we can access edge-neighboring active cells on same level to also
449 // access of an edge to the edges "children". These are found from looking
450 // at the corresponding edge of children of inactive edge neighbors.
451 for (unsigned int line_idx = 0; line_idx < n_raw_lines; ++line_idx)
452 {
453 if ((line_to_cells[line_idx].size() > 0) &&
454 line_to_inactive_cells[line_idx].size() > 0)
455 {
456 // We now have cells to add (active ones) and edges to which they
457 // should be added (inactive cells).
458 const Triangulation<3>::cell_iterator inactive_cell(
460 line_to_inactive_cells[line_idx][0][0],
461 line_to_inactive_cells[line_idx][0][1]);
462 const unsigned int neighbor_line =
463 line_to_inactive_cells[line_idx][0][2];
464
465 for (unsigned int c = 0; c < 2; ++c)
466 {
467 const auto &child =
468 inactive_cell->child(line_to_children[neighbor_line][c]);
469 const unsigned int child_line_idx =
470 child->line_index(neighbor_line);
471
472 Assert(child->is_active(), ExcInternalError());
473
474 // Now add all active cells
475 for (const auto &cl : line_to_cells[line_idx])
476 line_to_cells[child_line_idx].push_back(cl);
477 }
478 }
479 }
480 }
481
482
483
484 template <int dim>
485 inline std::vector<std::vector<bool>>
487 const ::hp::FECollection<dim> &fe_collection)
488 {
489 std::vector<std::vector<bool>> supported_components(
490 fe_collection.size(),
491 std::vector<bool>(fe_collection.n_components(), false));
492
493 for (unsigned int i = 0; i < fe_collection.size(); ++i)
494 {
495 for (unsigned int base_element_index = 0, comp = 0;
496 base_element_index < fe_collection[i].n_base_elements();
497 ++base_element_index)
498 for (unsigned int c = 0;
499 c < fe_collection[i].element_multiplicity(base_element_index);
500 ++c, ++comp)
501 if (dim == 1 ||
502 (dynamic_cast<const FE_Q<dim> *>(
503 &fe_collection[i].base_element(base_element_index)) ==
504 nullptr &&
505 dynamic_cast<const FE_Q_iso_Q1<dim> *>(
506 &fe_collection[i].base_element(base_element_index)) ==
507 nullptr))
508 supported_components[i][comp] = false;
509 else
510 supported_components[i][comp] = true;
511 }
512
513 return supported_components;
514 }
515
516
517
518 template <int dim>
519 template <typename CellIterator>
520 inline ConstraintKinds
522 const CellIterator &cell) const
523 {
524 // TODO: for simplex or mixed meshes: nothing to do
525 if ((dim == 3 && line_to_cells.empty()) ||
526 (cell->reference_cell().is_hyper_cube() == false))
528
529 if (cell->level() == 0)
531
532 const std::uint16_t subcell =
533 cell->parent()->child_iterator_to_index(cell);
534 const std::uint16_t subcell_x = (subcell >> 0) & 1;
535 const std::uint16_t subcell_y = (subcell >> 1) & 1;
536 const std::uint16_t subcell_z = (subcell >> 2) & 1;
537
538 std::uint16_t face = 0;
539 std::uint16_t edge = 0;
540
541 for (unsigned int direction = 0; direction < dim; ++direction)
542 {
543 const auto side = (subcell >> direction) & 1U;
544 const auto face_no = direction * 2 + side;
545
546 // ignore if at boundary
547 if (cell->at_boundary(face_no))
548 continue;
549
550 const auto &neighbor = cell->neighbor(face_no);
551
552 // ignore neighbors that are artificial or have the same level or
553 // have children
554 if (neighbor->has_children() || neighbor->is_artificial() ||
555 neighbor->level() == cell->level())
556 continue;
557
558 // Ignore if the neighbors are FE_Nothing
559 if (neighbor->get_fe().n_dofs_per_cell() == 0)
560 continue;
561
562 face |= 1 << direction;
563 }
564
565 if (dim == 3)
566 for (unsigned int direction = 0; direction < dim; ++direction)
567 if (face == 0 || face == (1 << direction))
568 {
569 const unsigned int line_no =
570 direction == 0 ?
571 (local_lines[0][subcell_y == 0][subcell_z == 0]) :
572 (direction == 1 ?
573 (local_lines[1][subcell_x == 0][subcell_z == 0]) :
574 (local_lines[2][subcell_x == 0][subcell_y == 0]));
575
576 const unsigned int line_index = cell->line_index(line_no);
577
578 const auto edge_neighbor =
579 std::find_if(line_to_cells[line_index].begin(),
580 line_to_cells[line_index].end(),
581 [&cell](const auto &edge_neighbor) {
583 &cell->get_triangulation(),
584 edge_neighbor[0],
585 edge_neighbor[1],
586 &cell->get_dof_handler());
587 return dof_cell.is_artificial() == false &&
588 dof_cell.level() < cell->level() &&
589 dof_cell.get_fe().n_dofs_per_cell() > 0;
590 });
591
592 if (edge_neighbor == line_to_cells[line_index].end())
593 continue;
594
595 edge |= 1 << direction;
596 }
597
598 if ((face == 0) && (edge == 0))
600
601 const std::uint16_t inverted_subcell = (subcell ^ (dim == 2 ? 3 : 7));
602
603 const auto refinement_configuration = static_cast<ConstraintKinds>(
604 inverted_subcell + (face << 3) + (edge << 6));
605 Assert(check(refinement_configuration, dim), ExcInternalError());
606 return refinement_configuration;
607 }
608
609
610
611 template <int dim>
612 template <typename CellIterator>
613 inline void
615 const CellIterator &cell,
616 const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner,
617 const std::vector<std::vector<unsigned int>> &lexicographic_mapping,
618 const std::vector<std::vector<bool>> &supported_components,
619 const ConstraintKinds &refinement_configuration,
620 std::vector<types::global_dof_index> &dof_indices) const
621 {
622 if (std::find(supported_components[cell->active_fe_index()].begin(),
623 supported_components[cell->active_fe_index()].end(),
624 true) ==
625 supported_components[cell->active_fe_index()].end())
626 return;
627
628 const auto &fe = cell->get_fe();
629 AssertDimension(fe.n_unique_faces(), 1);
630
631 std::vector<std::vector<unsigned int>>
632 component_to_system_index_face_array(fe.n_components());
633
634 for (unsigned int i = 0; i < fe.n_dofs_per_face(0); ++i)
635 component_to_system_index_face_array
636 [fe.face_system_to_component_index(i, /*face_no=*/0).first]
637 .push_back(i);
638
639 std::vector<unsigned int> idx_offset = {0};
640
641 for (unsigned int base_element_index = 0;
642 base_element_index < fe.n_base_elements();
643 ++base_element_index)
644 for (unsigned int c = 0;
645 c < fe.element_multiplicity(base_element_index);
646 ++c)
647 idx_offset.push_back(
648 idx_offset.back() +
649 fe.base_element(base_element_index).n_dofs_per_cell());
650
651 std::vector<types::global_dof_index> neighbor_dofs_all(idx_offset.back());
652 std::vector<types::global_dof_index> neighbor_dofs_all_temp(
653 idx_offset.back());
654 std::vector<types::global_dof_index> neighbor_dofs_face(
655 fe.n_dofs_per_face(/*face_no=*/0));
656
657
658 const auto get_face_idx = [](const auto n_dofs_1d,
659 const auto face_no,
660 const auto i,
661 const auto j) -> unsigned int {
662 const auto direction = face_no / 2;
663 const auto side = face_no % 2;
664 const auto offset = (side == 1) ? (n_dofs_1d - 1) : 0;
665
666 if (dim == 2)
667 return (direction == 0) ? (n_dofs_1d * i + offset) :
668 (n_dofs_1d * offset + i);
669 else if (dim == 3)
670 switch (direction)
671 {
672 case 0:
673 return n_dofs_1d * n_dofs_1d * i + n_dofs_1d * j + offset;
674 case 1:
675 return n_dofs_1d * n_dofs_1d * j + n_dofs_1d * offset + i;
676 case 2:
677 return n_dofs_1d * n_dofs_1d * offset + n_dofs_1d * i + j;
678 default:
680 }
681
683
684 return 0;
685 };
686
687 const std::uint16_t kind =
688 static_cast<std::uint16_t>(refinement_configuration);
689 const std::uint16_t subcell = (kind >> 0) & 7;
690 const std::uint16_t subcell_x = (subcell >> 0) & 1;
691 const std::uint16_t subcell_y = (subcell >> 1) & 1;
692 const std::uint16_t subcell_z = (subcell >> 2) & 1;
693 const std::uint16_t face = (kind >> 3) & 7;
694 const std::uint16_t edge = (kind >> 6) & 7;
695
696 for (unsigned int direction = 0; direction < dim; ++direction)
697 if ((face >> direction) & 1U)
698 {
699 const auto side = ((subcell >> direction) & 1U) == 0;
700 const auto face_no = direction * 2 + side;
701
702 // read DoFs of parent of face, ...
703 cell->neighbor(face_no)
704 ->face(cell->neighbor_face_no(face_no))
705 ->get_dof_indices(neighbor_dofs_face,
706 cell->neighbor(face_no)->active_fe_index());
707
708 // ... convert the global DoFs to serial ones, and ...
709 if (partitioner)
710 for (auto &index : neighbor_dofs_face)
711 index = partitioner->global_to_local(index);
712
713 for (unsigned int base_element_index = 0, comp = 0;
714 base_element_index < fe.n_base_elements();
715 ++base_element_index)
716 for (unsigned int c = 0;
717 c < fe.element_multiplicity(base_element_index);
718 ++c, ++comp)
719 {
720 if (supported_components[cell->active_fe_index()][comp] ==
721 false)
722 continue;
723
724 const unsigned int n_dofs_1d =
725 cell->get_fe()
726 .base_element(base_element_index)
727 .tensor_degree() +
728 1;
729 const unsigned int dofs_per_face =
730 Utilities::pow(n_dofs_1d, dim - 1);
731 std::vector<types::global_dof_index> neighbor_dofs(
732 dofs_per_face);
733 const auto lex_face_mapping =
735 n_dofs_1d - 1);
736
737 // ... extract the DoFs of the current component
738 for (unsigned int i = 0; i < dofs_per_face; ++i)
739 neighbor_dofs[i] = neighbor_dofs_face
740 [component_to_system_index_face_array[comp][i]];
741
742 // fix DoFs depending on orientation, rotation, and flip
743 if (dim == 2)
744 {
745 // TODO: this needs to be implemented for simplices but
746 // all-quad meshes are OK
747 Assert(cell->combined_face_orientation(face_no) ==
750 }
751 else if (dim == 3)
752 {
753 orient_face(cell->combined_face_orientation(face_no),
754 n_dofs_1d,
755 neighbor_dofs);
756 }
757 else
758 {
760 }
761
762 // update DoF map
763 for (unsigned int i = 0, k = 0; i < n_dofs_1d; ++i)
764 for (unsigned int j = 0; j < (dim == 2 ? 1 : n_dofs_1d);
765 ++j, ++k)
766 dof_indices[get_face_idx(n_dofs_1d, face_no, i, j) +
767 idx_offset[comp]] =
768 neighbor_dofs[lex_face_mapping[k]];
769 }
770 }
771
772 if (dim == 3)
773 for (unsigned int direction = 0; direction < dim; ++direction)
774 if ((edge >> direction) & 1U)
775 {
776 const unsigned int line_no =
777 direction == 0 ?
778 (local_lines[0][subcell_y][subcell_z]) :
779 (direction == 1 ? (local_lines[1][subcell_x][subcell_z]) :
780 (local_lines[2][subcell_x][subcell_y]));
781
782 const unsigned int line_index = cell->line(line_no)->index();
783
784 const auto edge_neighbor =
785 std::find_if(line_to_cells[line_index].begin(),
786 line_to_cells[line_index].end(),
787 [&cell](const auto &edge_array) {
789 edge_neighbor(&cell->get_triangulation(),
790 edge_array[0],
791 edge_array[1]);
792 return edge_neighbor->is_artificial() == false &&
793 edge_neighbor->level() < cell->level();
794 });
795
796 if (edge_neighbor == line_to_cells[line_index].end())
797 continue;
798
799 const DoFCellAccessor<dim, dim, false> neighbor_cell(
800 &cell->get_triangulation(),
801 (*edge_neighbor)[0],
802 (*edge_neighbor)[1],
803 &cell->get_dof_handler());
804 const auto local_line_neighbor = (*edge_neighbor)[2];
805
806 neighbor_cell.get_dof_indices(neighbor_dofs_all);
807
808 if (partitioner)
809 for (auto &index : neighbor_dofs_all)
810 index = partitioner->global_to_local(index);
811
812 for (unsigned int i = 0; i < neighbor_dofs_all_temp.size(); ++i)
813 neighbor_dofs_all_temp[i] = neighbor_dofs_all
814 [lexicographic_mapping[cell->active_fe_index()][i]];
815
816 const bool flipped =
817 cell->line_orientation(line_no) !=
818 neighbor_cell.line_orientation(local_line_neighbor);
819
820 for (unsigned int base_element_index = 0, comp = 0;
821 base_element_index < fe.n_base_elements();
822 ++base_element_index)
823 for (unsigned int c = 0;
824 c < fe.element_multiplicity(base_element_index);
825 ++c, ++comp)
826 {
827 if (supported_components[cell->active_fe_index()][comp] ==
828 false)
829 continue;
830
831 const unsigned int n_dofs_1d =
832 cell->get_fe()
833 .base_element(base_element_index)
834 .tensor_degree() +
835 1;
836
837 for (unsigned int i = 0; i < n_dofs_1d; ++i)
838 dof_indices[line_dof_idx(line_no, i, n_dofs_1d) +
839 idx_offset[comp]] = neighbor_dofs_all_temp
840 [line_dof_idx(local_line_neighbor,
841 flipped ? (n_dofs_1d - 1 - i) : i,
842 n_dofs_1d) +
843 idx_offset[comp]];
844 }
845 }
846 }
847
848
849
850 template <int dim>
851 template <typename CellIterator>
852 inline bool
854 const CellIterator &cell,
855 const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner,
856 const std::vector<std::vector<unsigned int>> &lexicographic_mapping,
857 std::vector<types::global_dof_index> &dof_indices,
858 const ArrayView<ConstraintKinds> &masks) const
859 {
860 // 1) check if finite elements support fast hanging-node algorithm
861 const auto supported_components = compute_supported_components(
862 cell->get_dof_handler().get_fe_collection());
863
864 if (std::none_of(supported_components.begin(),
865 supported_components.end(),
866 [](const auto &a) {
867 return *std::max_element(a.begin(), a.end());
868 }))
869 return false;
870
871 // 2) determine the refinement configuration of the cell
872 const auto refinement_configuration =
873 compute_refinement_configuration(cell);
874
875 if (refinement_configuration == ConstraintKinds::unconstrained)
876 return false;
877
878 // 3) update DoF indices of cell for specified components
879 update_dof_indices(cell,
880 partitioner,
881 lexicographic_mapping,
882 supported_components,
883 refinement_configuration,
884 dof_indices);
885
886 // 4) TODO: copy refinement configuration to all components
887 for (unsigned int c = 0; c < supported_components[0].size(); ++c)
888 if (supported_components[cell->active_fe_index()][c])
889 masks[c] = refinement_configuration;
890
891 return true;
892 }
893
894
895
896 template <int dim>
897 inline void
899 unsigned int &subface_index) const
900 {
901 const unsigned int rot_mapping[4] = {2, 0, 3, 1};
902
903 times = times % 4;
904 times = times < 0 ? times + 4 : times;
905 for (int t = 0; t < times; ++t)
906 subface_index = rot_mapping[subface_index];
907 }
908
909
910
911 template <int dim>
912 inline void
914 const types::geometric_orientation combined_orientation,
915 const unsigned int n_dofs_1d,
916 std::vector<types::global_dof_index> &dofs) const
917 {
918 const auto [orientation, rotation, flip] =
919 ::internal::split_face_orientation(combined_orientation);
920 const int n_rotations =
921 rotation || flip ? 4 - int(rotation) - 2 * int(flip) : 0;
922
923 const unsigned int rot_mapping[4] = {2, 0, 3, 1};
924 Assert(n_dofs_1d > 1, ExcInternalError());
925 // 'per line' has the same meaning here as in FiniteElementData, i.e., the
926 // number of dofs assigned to a line (not including vertices).
927 const unsigned int dofs_per_line = n_dofs_1d - 2;
928
929 // rotate:
930 std::vector<types::global_dof_index> copy(dofs.size());
931 for (int t = 0; t < n_rotations; ++t)
932 {
933 std::swap(copy, dofs);
934
935 // Vertices
936 for (unsigned int i = 0; i < 4; ++i)
937 dofs[rot_mapping[i]] = copy[i];
938
939 // Edges
940 unsigned int offset = 4;
941 for (unsigned int i = 0; i < dofs_per_line; ++i)
942 {
943 // Left edge
944 dofs[offset + i] =
945 copy[offset + 2 * dofs_per_line + (dofs_per_line - 1 - i)];
946 // Right edge
947 dofs[offset + dofs_per_line + i] =
948 copy[offset + 3 * dofs_per_line + (dofs_per_line - 1 - i)];
949 // Bottom edge
950 dofs[offset + 2 * dofs_per_line + i] =
951 copy[offset + dofs_per_line + i];
952 // Top edge
953 dofs[offset + 3 * dofs_per_line + i] = copy[offset + i];
954 }
955
956 // Interior points
957 offset += 4 * dofs_per_line;
958
959 for (unsigned int i = 0; i < dofs_per_line; ++i)
960 for (unsigned int j = 0; j < dofs_per_line; ++j)
961 dofs[offset + i * dofs_per_line + j] =
962 copy[offset + j * dofs_per_line + (dofs_per_line - 1 - i)];
963 }
964
965 // transpose (note that we are using the standard geometric orientation
966 // here so orientation = true is the default):
967 if (!orientation)
968 {
969 copy = dofs;
970
971 // Vertices
972 dofs[1] = copy[2];
973 dofs[2] = copy[1];
974
975 // Edges
976 unsigned int offset = 4;
977 for (unsigned int i = 0; i < dofs_per_line; ++i)
978 {
979 // Right edge
980 dofs[offset + i] = copy[offset + 2 * dofs_per_line + i];
981 // Left edge
982 dofs[offset + dofs_per_line + i] =
983 copy[offset + 3 * dofs_per_line + i];
984 // Bottom edge
985 dofs[offset + 2 * dofs_per_line + i] = copy[offset + i];
986 // Top edge
987 dofs[offset + 3 * dofs_per_line + i] =
988 copy[offset + dofs_per_line + i];
989 }
990
991 // Interior
992 offset += 4 * dofs_per_line;
993 for (unsigned int i = 0; i < dofs_per_line; ++i)
994 for (unsigned int j = 0; j < dofs_per_line; ++j)
995 dofs[offset + i * dofs_per_line + j] =
996 copy[offset + j * dofs_per_line + i];
997 }
998 }
999
1000
1001
1002 template <int dim>
1003 inline unsigned int
1005 unsigned int dof,
1006 unsigned int n_dofs_1d) const
1007 {
1008 unsigned int x, y, z;
1009
1010 const unsigned int fe_degree = n_dofs_1d - 1;
1011
1012 if (local_line < 8)
1013 {
1014 x = (local_line % 4 == 0) ? 0 :
1015 (local_line % 4 == 1) ? fe_degree :
1016 dof;
1017 y = (local_line % 4 == 2) ? 0 :
1018 (local_line % 4 == 3) ? fe_degree :
1019 dof;
1020 z = (local_line / 4) * fe_degree;
1021 }
1022 else
1023 {
1024 x = ((local_line - 8) % 2) * fe_degree;
1025 y = ((local_line - 8) / 2) * fe_degree;
1026 z = dof;
1027 }
1028
1029 return n_dofs_1d * n_dofs_1d * z + n_dofs_1d * y + x;
1030 }
1031
1032
1033
1034 template <int dim>
1035 void
1037 {
1038 if (subface == 1)
1039 subface = 2;
1040 else if (subface == 2)
1041 subface = 1;
1042 }
1043 } // namespace MatrixFreeFunctions
1044} // namespace internal
1045
1046
1048
1049#endif
bool is_active() const
const FiniteElement< dimension_, space_dimension_ > & get_fe() const
void get_dof_indices(std::vector< types::global_dof_index > &dof_indices) const
Definition fe_q.h:554
int index() const
int level() const
unsigned int line_index(const unsigned int i) const
ConstraintKinds compute_refinement_configuration(const CellIterator &cell) const
void setup_line_to_cell(const Triangulation< dim > &triangulation)
void orient_face(const types::geometric_orientation combined_orientation, const unsigned int n_dofs_1d, std::vector< types::global_dof_index > &dofs) const
std::vector< boost::container::small_vector< std::array< unsigned int, 3 >, 6 > > line_to_cells
static std::vector< std::vector< bool > > compute_supported_components(const ::hp::FECollection< dim > &fe)
HangingNodes(const Triangulation< dim > &triangualtion)
bool setup_constraints(const CellIterator &cell, const std::shared_ptr< const Utilities::MPI::Partitioner > &partitioner, const std::vector< std::vector< unsigned int > > &lexicographic_mapping, std::vector< types::global_dof_index > &dof_indices, const ArrayView< ConstraintKinds > &mask) const
void transpose_subface_index(unsigned int &subface) const
unsigned int line_dof_idx(int local_line, unsigned int dof, unsigned int n_dofs_1d) const
void update_dof_indices(const CellIterator &cell, const std::shared_ptr< const Utilities::MPI::Partitioner > &partitioner, const std::vector< std::vector< unsigned int > > &lexicographic_mapping, const std::vector< std::vector< bool > > &component_mask, const ConstraintKinds &refinement_configuration, std::vector< types::global_dof_index > &dof_indices) const
const ::ndarray< unsigned int, 3, 2, 2 > local_lines
void rotate_subface_index(int times, unsigned int &subface_index) const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:513
#define DEAL_II_HOST_DEVICE
Definition config.h:114
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:514
unsigned int cell_index
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInternalError()
#define DEAL_II_NOT_IMPLEMENTED()
std::size_t size
Definition mpi.cc:734
std::vector< unsigned int > lexicographic_to_hierarchic_numbering(unsigned int degree)
constexpr T pow(const T base, const int iexp)
Definition utilities.h:966
bool operator<(const ConstraintKinds f1, const ConstraintKinds f2)
ConstraintKinds decompress(const compressed_constraint_kind kind_in, const unsigned int dim)
bool check(const ConstraintKinds kind_in, const unsigned int dim)
compressed_constraint_kind compress(const ConstraintKinds kind_in, const unsigned int dim)
ConstraintKinds & operator|=(ConstraintKinds &f1, const ConstraintKinds f2)
std::size_t memory_consumption(const ConstraintKinds &)
std::uint8_t compressed_constraint_kind
Definition dof_info.h:86
ConstraintKinds operator&(const ConstraintKinds f1, const ConstraintKinds f2)
bool operator!=(const ConstraintKinds f1, const ConstraintKinds f2)
constexpr compressed_constraint_kind unconstrained_compressed_constraint_kind
ConstraintKinds operator|(const ConstraintKinds f1, const ConstraintKinds f2)
std::tuple< bool, bool, bool > split_face_orientation(const types::geometric_orientation combined_face_orientation)
constexpr types::geometric_orientation default_geometric_orientation
Definition types.h:346
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
unsigned char geometric_orientation
Definition types.h:40
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation