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